# Workshop Programme

## for period 24 - 28 September 2012

### Solution of Partial Differential Equations on the Sphere (PDEs on the sphere)

24 - 28 September 2012

Timetable

Monday 24 September | ||||

09:00-10:20 | Registration | |||

09:50-10:20 | Morning Coffee | |||

10:20-10:30 | Welcome from the Deputy Director | |||

Session: Horizontal Session | ||||

10:30-10:55 | John Thuburn (University of Exeter) |
|||

Mimetic Semi-implicit Solution of the Shallow Water Equations on Hexagonal-Icosahedral and Cubed-Sphere Grids | Sem 1 | |||

A new algorithm is presented for the solution of the shallow water equations on quasi-uniform spherical grids. It combines a mimetic spatial discretization with a Crank-Nicolson time scheme for fast waves and an accurate and conservative forward-in-time advection scheme for mass and potential vorticity. The algorithm is tested on two families of grids: hexagonal-icosahedral Voronoi grids, and modified equiangular cube-sphere grids. For the cubed-sphere case, a key ingredient is the development of a suitable discrete Hodge star operator for the non-orthogonal grid. Results of several test cases will be presented. The results confirm a number of desirable properties for which the scheme was designed: exact mass conservation, very good available energy and potential enstrophy conservation, vanishing curl of grad, steady geostrophic modes, and accurate PV advection. The scheme is stable for large wave Courant numbers and for advective Courant numbers up to about one. The accuracy of the scheme appears to be limited by the accuracy of the various mimetic spatial operators. On the hexagonal grid there is no evidence for damaging effects of computational Rossby modes, despite attempts to force them explicitly. |
||||

10:55-11:20 | Thomas Dubos (École Polytechnique) |
|||

Variational derivation of energy-conserving finite-difference schemes for geophysical fluid equations | Sem 1 | |||

At the continuous level, the so-called vector-invariant form of the equations of fluid motion can be obtained from Hamilton's principle of least action, using the Euler-Poincaré formalism. In this form of the equations, the mass flux and Bernoulli function appear as functional derivatives of the Hamiltonian. The integral conservation of energy and the Lagrangian conservation of potential vorticity then follow straightforwardly. This setting can be imitated at the discrete level to yield energy-conserving schemes. Key ingredients include an energy-conserving vector product, a discrete rule of integration by parts, a discrete approximation of total energy, and the definition of the discrete mass flux and Bernoulli function from partial derivatives of the discrete energy. For the rotating shallow-water equations, it turns out that schemes by Sadourny (1975) on Cartesian meshes and by Bonaventura & Ringler (2005) on Delaunay-Voronoi meshes fit in the above framework. Furthermore new schemes can be obtained. For instance it is possible to modify the discrete mass flux to match a modification of the discrete energy, as Renner (1981) and Skamarock et al. (2012) did to suppress a numerical instability. Finally, a so-called discrete Hodge star operator is obtained on a non-orthogonal pair of primal and dual meshes. This operator is part of a potential-vorticity conserving scheme of the shallow-water equations on meshes for which an orthogonal dual does not exist, like the equiangular cubed sphere. Potential application to sets of three-dimensional equations will be discussed. |
||||

11:20-11:45 | Colin Cotter (Imperial College London) |
|||

Finite element exterior calculus framework for geophysical fluid dynamics | Sem 1 | |||

We present a finite element framework that extends the mimetic approach of Ringler, Thuburn, Skamarock and Klemp (2010) to finite element methods. The framework preserves the following benefits of the mimetic approach: local mass conservation, steady geostrophic linear modes on the f-plane, consistent conservative potential vorticity advection, but allows for selection of higher-order finite element discretisations and adjustment of the balance of vorticity and mass degrees of freedom. We will explain the foundations of the framework and present some preliminary numerical results which form part of the UK Met Office/NERC/STFC GungHo dynamical core project. |
||||

11:45-12:10 | Stefan Vater (Freie Universität Berlin) |
|||

A multilevel time integrator for computing longwave shallow water flows at low Froude numbers | Sem 1 | |||

A new multilevel semi-implicit scheme for the computation of low Froude number shallow water flows is presented. Motivated by the needs of atmospheric flow applications, it aims to minimize dispersion and amplitude errors in the computation of long wave gravity waves. While it correctly balances "slaved" dynamics of short-wave solution components induced by slow forcing, the method eliminates freely propagating compressible short-wave modes, which are under- resolved in time. This is achieved through the multilevel approach borrowing ideas from multigrid schemes for elliptic equations. The scheme is second order accurate and admits time steps depending only on the flow velocity. It incorporates a predictor step using a Godunov-type method for hyperbolic conservation laws and two elliptic corrections. The multilevel method is initially derived for the one-dimensional linearized shallow water equations. Scale-wise decomposition of the data enables a scale- dependent blending of time integrators with different principal features. To guide the selection of these integrators, the discrete-dispersion relations of some standard second-order schemes are analyzed, and their response to high wave number low frequency source terms is discussed. The resulting method essentially consists of the solution of a Helmholtz problem on the original fine grid, where the operator and the right hand side incorporate the multiscale information of the discretization. The performance of the method in the linear case is illustrated on a test case with "multiscale" initial data and a problem with a slowly varying high wave number source term. For the computation of fully nonlinear shallow water flows, a projection method for zero Froude number flows is generalized by incorporating the local time derivatives of the height. This semi-implicit method is combined with the multilevel method for the linearized equations to obtain the above mentioned properties of the scheme. Numerical tests address the scheme's ability to correctly cover the asymptotic flow regime of long-wave gravity waves passing over short-range topography and its balancing properties for the lake at rest. |
||||

14:00-14:25 | Sergey Danilov (Alfred-Wegener-Institut für Polar- und Meeresforschung (AWI)) |
|||

Multiscale ocean simulations with FESOM | Sem 1 | |||

Finite-Element Sea-ice Ocean circulation Model (FESOM) uses unstructured triangular meshes and is employed either as standing-alone global ocean model or as part of FESOM-ECHAM5/6 coupled setup. Several examples pertaining to studies focused on polar oceans will be given to characterize the potential of FESOM and also to indicate typical difficulties related to the strongly variable resolution. Brief intercomparison of FESOM to a global cell-vertex finite-volume setup is presented. While it is more numerically efficient than the finite-element setup of FESOM, it requires special measures to maintain stable performance because of its too large velocity space. Most effective stabilization is provided by either computing the momentum advection on scalar (median-dual) control volumes, or making use of the vector-invariant form. Both approaches effectively trim the nonlinear term in the momentum equation to the size of scalar space. |
||||

14:25-14:50 | Jürgen Steppeler (Deutscher Wetterdienst (DWD)) |
|||

Uniformly third Order conserving Schems on Polygonal Grids | Sem 1 | |||

Uniformly third Order conserving Schems on Polygonal Grids. The interest in polygonal grids is increasing. They are an alternative to the more commonly used spectral and latitude longitude grids. Among other advantages they offer the possibility of a rather uniform cover of the sphere with grid cells. Other advantages concern the ease of using multiprocessing computers and using special vertical treatments, such as shaved cells. Well known examples of polygonal grids are the cube sphere and the icosahedral grid. After initial research by Sadourny and Williamsson the practicability of this approach was shown by Baumgardner and Steppeler. In particular Baumgardner showed that problems with some approaches can be traced back to the fact that for slightly irregular resolution methods are not uniformly second order. After correcting this problem Baumgardner was able to show that problems arising from irregular grids do not occur. Steppeler generalized this approach to third order. Both Baumgardners and Steppelers approaches were non conser ving. A generalization to conserving schemes will be presented and computational examples given. Another high order approach is the pecral element method, which currently is available for orders 4 an higher only. The approach presented can be considered as a version of third order spectral elements. The advantages of third order schemes over even higher order approaches will be discussed. |
||||

14:50-15:15 | Ramachandran Nair (National Center for Atmospheric Research) |
|||

A Semi-Lagrangian Discontinuous Galerkin (SLDG) Conservative Transport Scheme on the Cubed-Sphere | Sem 1 | |||

The discontinuous Galerkin (DG) method combines fine features of high-order accurate finite-element and finite-volume methods. Because of its geometric flexibility and high parallel efficiency, DG method is becoming increasingly popular in atmospheric and ocean modeling. However, a major drawback of DG method is its stringent CFL stability restriction associated with explicit time-stepping. A way to get around this issue is to combine DG method with a Lagrangian approach based on the characteristic Lagrange-Galerkin philosophy. Unfortunately, a fully 2D approach combining DG and Lagrangian methods is algorithmically complex and computationally expensive for practical application, particularly for non-orthogonal curvilinear geometry such as the cubed-sphere grid system. We adopt a dimension-splitting approach where a regular semi-Lagrangian (SL) scheme is combined with the DG method. The resulting SLDG scheme employs a sequence of 1D operations for solving transport equation on the cubed-sphere. The SLDG scheme is inherently conservative and has the option to incorporate a local positivity-preserving filter for tracers. A novel feature of the SLDG algorithm is that it can be used for multi-tracer transport for global models employing spectral-element (structured or unstructured) grids. The SLDG scheme is tested for various benchmark advection test-suites on the sphere and results will be presented in the seminar. |
||||

15:45-16:10 | Jean-Pierre Croisille (Université de Lorraine) |
|||

Compact finite difference schemes on the Cubed-Sphere | Sem 1 | |||

The Cubed-Sphere is a spherical grid made of six quasi-cartesian square like patches. It was originally introduced by Sadourny some forty years ago. We extend to this grid the design of high-order finite difference compact operators. Such discrete operators are used in Computational Fluid Dynamics on structured grids for applications such as Direct Numerical Simulation of turbulent flows, or aeroacoustics problems. We consider in this work the design of a uniformly fourth-order accurate spherical gradient. The main approximation principle consists in defining a network of great circles covering the Cubed-Sphere along which a high-order hermitian gradient can be calculated. This procedure allows a natural treatment at the interface of the six patches. The main interest of the approach is a fully symmetric approximation system on the Cubed-Sphere. We numerically demonstrate the accuracy of the approximate gradient on several test problems, in particular the cosine-bell test-case of Williamson et al. for climatology. |
||||

16:10-16:35 | Nicholas Kevlahan (McMaster University) |
|||

A conservative adaptive wavelet method for the shallow water equations on staggered grids | Sem 1 | |||

This paper presents the first dynamically adaptive wavelet method for the shallow water equations on a staggered hexagonal C-grid. Pressure is located at the centres of the primal grid (hexagons) and velocity is located at the edges of the dual grid (triangles). Distinct biorthogonal second generation wavelet transforms are developed for the pressure and the velocity. These wavelet transforms are based on second-order accurate interpolation and restriction operators. Together with compatible restriction operators for the mass flux and Bernoulli function, they ensure that mass is conserved and that there is no numerical generation of vorticity when solving the shallow water equations. Grid refinement relies on appropriate thresholding of the wavelet coefficients, allowing error control in both the quasi-geostrophic and inertia-gravity wave regimes. The shallow water equations are discretized on the dynamically adapted multiscale grid using a mass and potential-enstrophy conserving finite-difference scheme. The conservation and error control properties of the method are verified by applying it to a propagating inertia-gravity wave packet and to rotating shallow water turbulence. Significant savings in the number of degrees of freedom are achieved even in the case of rotating shallow-water turbulence. The numerical dissipation introduced by the grid adaptation is quantified. The method has been designed so it can be extended easily to the icosahedral subdivision of the sphere. This work provides important building blocks for the development of fully adaptive general circulation models. |
||||

16:35-17:00 | Tae-Jin Oh (Korea Institute of Atmospheric Prediction Systems (KIAPS)) |
|||

Dynamical Core Developments at KIAPS | Sem 1 | |||

Korea Institute of Atmospheric Prediction Systems (KIAPS) is a new organization founded to develop the next generation operational numerical weather prediction (NWP) model for Korea Meteorological Administration (KMA). With the increasing demand for global high-resolution simulation, scalability becomes an important issue and highly scalable numerical methods such as spectral element (continuous Galerkin, CG) or discontinuous Galerkin (DG) methods are gaining interest. Although conventional finite difference (FD) or finite volume (FV) methods offer excellent efficiency, it is difficult to formulate high order schemes on nonorthogonal grid structures which is needed to avoid grid singularities. On the other hand, CG/DG methods are not constrained to lower order on unstructured, nonorthogonal grids. We present high order convergence properties of CG/DG methods for advection, shallow water equations on structured and/or unstructured grids in one and/or two dimensions. In order t o match the high order spatial truncation error, an explicit general order m-stage m–1 order strong stability preserving (SSP) Runge-Kutta time integrator is used. With this setup, we can obtain arbitrary order convergence rates. |
||||

12:30-13:30 | Lunch at Wolfson Court | |||

13:30-14:00 | Poster Session | |||

15:15-15:45 | Afternoon Tea | |||

17:00-18:00 | Welcome Drinks Reception and Poster Session |

Tuesday 25 September | ||||

Session: Advection Session | ||||

09:00-09:25 | Sarvesh Kumar Dubey (Indian Institute of Technology) |
|||

Slope-limited transport schemes using icosahedral hexagonal grid | Sem 1 | |||

In this work two simple advection schemes for unstructured meshes are proposed and implemented over spherical icosahedral-hexagonal grids. One of them is fully discrete in space and time while the other one is a semi discrete scheme with third order RungeKutta time integration. Both schemes use cell-wise linear reconstruction. We therefore also present two possible candidates for consistent gradient discretization over general grids. These gradients are designed in a finite volume sense with an adequate modification to guarantee convergence in the absence of a special grid optimization. Monotonicity of the advection schemes is enforced by a slope limiter, at contrast with the widely used of posterior approach of flux-corrected transport (FCT). Convergence of the proposed gradient reconstruction operators is verified numerically. It is found that the proposed modification is indeed necessary for convergence on non-optimized grids. Recently proposed advection test cases are used to evaluate the performance of the slope-limited advection schemes. It is verified that they are convergent and positive. We also compare these schemes to a variant where positivity is enforced by the FCT approach. FCT produces slightly less diffusion but it seems to be at the price of some nonphysical anti-diffusion. These results suggest that the proposed slope limited advection schemes are a viable option for icosahedral-hexagonal grids over sphere. |
||||

09:25-09:50 | Oswald Knoth (Leibniz Institute for Tropospheric Research, Leipzig) |
|||

Numerical Solution of the Advection Equation on Unstructured Spherical Grids with Logarithmic Reconstruction | Sem 1 | |||

There are numerous approaches for solving hyperbolic differential equations in the context of finite volume methods. One popular approach is the limiter free Local-Double-Logarithmic-Reconstruction (LDLR) of Artebrant and Schroll. The aim of this work is to construct a three-dimensional reconstructing function based on the LDLR for solving the advection equation on unstructured spherical grids. The new method should preserve the characteristics of the LDLR. That means in particular a reconstruction without use of limiters and with a small stencil of only the nearest neighbors of a particular cell. Also local extrema should be conserved while the local variation of the reconstruction within one cell should be under control. We propose an ansatz which works on unstructured polyhedral grids. To come up to this, an ansatz function with one logarithmic expression for each face of the polyghedron is constructed. Required gradients at cell face midpoints are determined by use of the Multi-Point-Flux-Approximation (MPFA) method. Further derivative information are obtained with the help of special barycentric coordinates. All necessary integrals of the ansatz functions can be computed exactly. The spatially discretized equations are combined with explicit Runge-Kutta methods to advance the solution in time. The new advection procedure is numerically evaluated with standard test cases from the literature on different unstructured spherical grids. |
||||

09:50-10:15 | Janakiraman Subburathnam (Centre for Development of Advanced Computing, (C-DAC)) |
|||

Linear advection characteristics of a variable resolution global spectral method on the sphere | Sem 1 | |||

Advection experiments are conducted with a variable resolution global spectral method on the sphere. The variable resolution grid has finer resolution over the tropics and the resolution decreases smoothly as we move towards the poles [Janakiraman et al (2012) ]. An Eulerian formulation of the linear advection is used to for the spectral discretization. Near dispersion-free advection is achieved on the high resolution tropical belt. The transport across homogeneous resolution regions produce very less dispersion errors. Transport over the poles result in severe grid representation errors. It is shown that increasing the resolution over the reatly reduces this error. Transport of a feature from apoint close to poles but not over it does not produce such representation errors. A comparison of time-schemes such as 4-th order Runge-Kutta method, Magazenkov scheme and Leap-frog scheme for the advection experiment is also presented. Reference(s) : S. Janakiraman, Ravi S Nanjundiah and A.S. Vasudeva Murthy, A novel variable resolution global spectral method on the sphere, Vol. 231, No. 7, pp 2794--2810, 2012. |
||||

10:15-10:40 | Kara Peterson (Sandia National Laboratories) |
|||

Optimization-based Conservative Transport on the Sphere | Sem 1 | |||

We present a new optimization-based conservative transport algorithm for scalar quantities (i.e. mass) that preserves monotonicity without the use of flux limiters. The method is formulated as a singly linearly constrained quadratic program with simple bounds where the net mass updates to the cell are the optimization variables. The objective is to minimize the discrepancy between a target or high-order mass update and a mass update that satisfies physical bounds. In this way, we separate accuracy considerations, handled by the objective functional from the enforcement of physical bounds, handled by the contraints. With this structure mass conservation is incorporated as a constraint and a simple, efficient, and easily parallelizable optimization algorithm is obtained. This algorithm has been extended to a latitude-longitude grid for two-dimensional remapping and transport on the sphere. Results for several standard test problems on the sphere will be shown to illustrate the accuracy and robustness of the method. |
||||

10:40-11:10 | Morning Coffee | |||

Session: Adaptive Meshing | ||||

11:10-11:35 | Colin Zarzycki (University of Michigan) |
|||

Improving tropical cyclone representation in general circulation models through the use of variable resolution | Sem 1 | |||

Modeling of tropical cyclones in General Circulation Models (GCMs) has been a historically difficult task due to issues such as relatively small storm sizes and intense convective processes. However, recent advances in GCM model design coupled with improvements in computing ability now allow for GCM simulations with grid spacings as small as 15-30 km. This presentation evaluates the potential of GCMs at these high horizontal resolutions to simulate tropical cyclones. In particular, we explore a novel variable-resolution mesh approach that allows for high spatial resolutions in areas of interest, such as low-latitude ocean basins where tropical cyclones are prevalent. Such GCM designs with variable-resolution meshes have the potential to become a future tool in weather forecasting as well as for regional climate assessments. A statically-nested, variable-resolution mesh option has recently been introduced into the National Center for Atmospheric Research (NCAR) Community Atmosphere Model's (CAM) Spectral Element (SE) dynamical core. We present preliminary CAM-SE model simulations using an idealized tropical cyclone test case with a variety of grid sizes and refinement scales. We evaluate the evolution of tropical cyclones initialized at various locations in or near grid scale transition regions. Specific focus is centered on factors crucial to storm cyclogenesis and maintenance such as air-sea interaction and vertical development. In addition to short-term, deterministic tests, we also investigate the performance of multi-resolution meshes in longer-term climate simulations, including attention paid to the dependance of non-seeded tropical cyclone genesis on spatial resolution and the subsequent implications for regional climate modeling within a global modeling framework. We also discuss pot ential computational consequences of using such a setup in either process or climate studies. |
||||

11:35-12:00 | Phillip Colella (Lawrence Berkeley National Laboratory) |
|||

Adaptive High-order Finite Volume Discretizations on Spherical Thin Shells | Sem 1 | |||

We present an adaptive, conservative finite volume approach applicable to solving hyperbolic PDE's on both 2D surface and 3D thin shells on the sphere. The starting point for this method is the equiangular cubed-sphere mapping, which maps six rectangular coordinate patches (blocks) onto the sphere. The images of these blocks form a disjoint union covering the sphere, with the mappings of adjacent blocks being continuous, but not differentiable, at block boundaries. Our method uses a fourth-order accurate discretization to compute flux averages on faces, with a higher-order least squares-based interpolation to compute stencil operations near block boundaries. To suppress oscillations at discontinuities and underresolved gradients, we use a limiter that preserves fourth-order accuracy at smooth extrema, and a redistribution scheme to preserve positivity where appropriate for advected scalars. By using both space- and time-adaptive mesh refinement, the solver allocates comp utational effort only where greater accuracy is needed. The resulting method is demonstrated to be fourth-order accurate for advection and shallow water equation model problems, and robust at solution discontinuities. We will also present an approach for the compressible Euler equations on a 3D thin spherical shell. Refinement is performed only in the horizontal directions, The radial direction is treated implicitly (using a fourth-order RK IMEX scheme) to eliminate time step constraints from vertical acoustic waves. |
||||

12:00-12:25 | Jonathan Lambrechts (Université Catholique de Louvain) |
|||

Generation of Provably Correct Curvilinear Meshes | Sem 1 | |||

The development of high-order numerical technologies for CFD is underway for many years now. For example, Discontinuous Galerkin methods (DGM) have been largely studied in the literature, initially in a quite theoretical context, and now in the application point of view. In many contributions, it is shown that the accuracy of the method strongly depends of the accuracy of the geometrical discretization. In other words, the following question is raised:we have the high order methods, but how do we get the meshes? This talk focus on the generation of highly curved ocean meshes of polynomial order 2 to 4. In the first part, we propose a robust procedure that allows to build a curvilinear mesh for which every element is guaranteed to be valid. The technique builds on standard optimization method (BICG) combined with a log-barrier objective function to guarantee the positivity of the elements Jacobian and thus the validity of the elements. To be valid is not the only requirement for a good-quality mesh. If the temporal discretization is explicit, even a valid element can lead to a very stringent constraint on the stable time step. The second part of the talk is devoted to the optimization of the curvilinear ocean meshes to obtain large stable time steps. |
||||

14:00-14:25 | Donna Calhoun (Boise State University) |
|||

A logically Cartesian, adaptively refined two-patch sphere grid for modeling transport in the atmosphere | Sem 1 | |||

Recently, we demonstrated results using a second order finite volume scheme on a novel, logically Cartesian, two-patch 2d sphere grid. The numerical scheme, based on the wave-propagation algorithms in Clawpack ( R. J. LeVeque, Univ. of Washington), was easily adapted to the single grid Cartesian layout of our two-patch sphere grid mapping. Furthermore, we were able to use an existing adaptive mesh refinement patch-based (AMR) code (Berger, Oliger et al.) to run computational efficient simulations. In our current work, we are developing a new AMR code which uses wave propagation algorithms on non-overlapping AMR grids stored as leaves in a forest of quad- or oct-trees. The underlying tree-based code, p4est (Carsten Burstedde, Univ. of Bonn) manages the multi-block connectivity in a multi-processor environment and has been shown to be highly scalable in applications of interest to geophysicists. This new AMR code, which we call ForestClaw, will easily handle the adaptivity for our two-patch sphere grid, as well as the cubed-sphere, and more generally, any multi-block geometry. We will present preliminary results from our efforts to use ForestClaw for modeling volcanic ash dispersal in the atmosphere. This is joint work with Carsten Burstedde (Univ. of Bonn) and researchers at the Cascade Volcanic Observatory (Vancouver, Washington). |
||||

14:25-14:50 | Andreas Müller (Naval Postgraduate School) |
|||

Accuracy of adaptive discontinuous Galerkin simulations | Sem 1 | |||

Adaptive mesh refinement generally serves to increase computational efficiency without compromising the accuracy of the numerical solution significantly. However it is an open question in which regions the spatial resolution can actually be coarsened without affecting the accuracy of the result significantly. This question is investigated for a specific example of dry atmospheric convection, namely the simulation of warm air bubbles. For this purpose a novel numerical model is developed that is tailored towards this specific meteorological problem. The compressible Euler equations are solved with a Discontinuous Galerkin method. Time integration is done with a semi-implicit approach and the dynamic grid adaptivity uses space filling curves via the AMATOS function library. The numerical model is validated with a convergence study and five standard test cases. A method is introduced which allows one to compare the accuracy between different choices of refinement regions even in a case when the exact solution is not known. Essentially this is done by comparing features of the solution that are strongly sensitive to spatial resolution. For a rising warm air bubble the additional error by using adaptivity is smaller than 1% of the total numerical error if the average number of elements used for the adaptive simulation is about 50% smaller than the number used for the simulation with the uniform fine-resolution grid. Correspondingly the adaptive simulation is almost two times faster than the uniform simulation. |
||||

14:50-15:15 | Giovanni Tumolo (Abdus Salam International Centre for Theoretical Physics) |
|||

A semi-implicit, semi-Lagrangian, p-adaptive Discontinuous Galerkin method for the rotating shallow water equations | Sem 1 | |||

As a first step towards construction and analysis of a DG based dynamical core for high resolution atmospheric modeling, a semi-implicit and semi-Lagrangian Discontinuous Galerkin method for the shallow water equations with rotation is proposed and analyzed. The method is equipped with a simple p-adaptivity criterion, that allows to adjust dynamically the number of local degrees of freedom employed to the local structure of the solution. Numerical results in the framework of two dimensional test cases prove that the method captures accurately and effectively the main features of linear gravity and inertial gravity waves. Also the solution of nonlinear Stommel problem is correctly simulated. The effectiveness of the method is also demonstrated by numerical results obtained at high Courant numbers and with automatic choice of the local approximation degree. The present research has been carried out during the PhD of the autor (supervisors F. Giorgi and L. Bonaventura) with financial support from the {\it Abdus Salam International Center for Theoretical Physics} and in collaboration with L.Bonaventura and M. Restelli {\it MOX - Politecnico di Milano}. See MOX-Report 04/2012 Tumolo, G.; Bonaventura, L.; Restelli, M. A semi-implicit, semi-Lagrangian, p-adaptive Discontinuous Galerkin method for the shallow water equations. |
||||

15:45-16:10 | Peter Bosler (University of Michigan) |
|||

Particle methods for geophysical flow on the sphere | Sem 1 | |||

We develop a fluid dynamics solver for spherical domains based on the Lagrangian form of the equations of motion and a tree-structured mesh of fluid particles. Initial discretizations are based on the cubed sphere or icosahedral triangles, but these arrangements distort as the particles move with the fluid velocity. A remeshing scheme based on interpolation of the Lagrangian parameter is applied at regular intervals to maintain computational accuracy as the flow evolves. Ongoing work with adaptive mesh refinement is introduced. We present solutions of the advection equation with prescribed velocities from recent test cases (Nair and Lauritzen, 2010), and solutions of the barotropic vorticity equation where velocity is given by an approximate Biot-Savart integral and midpoint rule quadrature for test cases including Rossby-Haurwitz waves and Gaussian vortices . The discrete equations in the barotropic vorticity equation are those of point vortices on the sphere (e.g. Bogo molov, 1977). For the shallow water equations, we must also include the effects of point sources upon the fluid particles; we discuss the challenges posed by divergent flows in this Lagrangian context and strategies for evaluating nonlinear forcing terms on irregular meshes. |
||||

16:10-16:35 | John Boyd (University of Michigan) |
|||

Progress in Radial Basis Function Methods: Adaptive Vortex-RBF Methods for the Sphere and Other Advances | Sem 1 | |||

Vortex methods, in which flow is modelled by overlapping Lagrangian vortex blobs, and radial basis function pseudospectral methods, which have proven their worth for irregular, adaptive grids, are combined into a single algorithm to compute flows on a sphere. The strengths and limitations of vortex-RBF methods for vortex-dominated flows with and without rotation will be illustrated by case studies including finite amplitude Rossby waves, vortex merger and deformation, and roll-up and instability. Other pertinent advances will be described as time permits. |
||||

16:35-17:00 | Michal Kopera (Naval Postgraduate School) |
|||

Adaptive mesh refinement for discontinuous Galerkin method on quadrilateral non-conforming grid | Sem 1 | |||

In recent years there has been a significant interest in using the discontinuous Galerkin method (DG) for solving fluid dynamics problems. Studies in [1], [2], and [3] have shown that the DG method is indeed a good choice for the construction of future non-hydrostatic numerical weather prediction models. It combines high-order accuracy of the solution with geometric flexibility of unstructured grids and exhibits excellent scaling properties. In order to increase the scale resolution capabilities of DG, as well as to take better advantage of available computing power, the use of adaptive mesh refinement (AMR) for a quadrilateral non-conforming grid is being investigated. The results of AMR implementation to a 2D version of the Non-hydrostatic Unified Model of the Atmosphere (NUMA) will be presented during the talk. Both static and dynamic h-adaptivity (element grid adaptivity) will be considered. Since increased local grid resolution has to impose severe constraints on the time-step of an explicit method, implicit-explicit (IMEX) and multi rate time integration will be discussed as well. As the latest trend in scientific computing is to employ clusters of GPUs for high-performance computations, the implementation of CUDA kernels into AMR NUMA will be examined. [1] Giraldo, F. & Restelli, M. (2008). A study of spectral element and discontinuous Galerkin methods for the Navier-Stokes equations in non-hydrostatic mesoscale atmospheric modelling: Equation sets and test cases. Journal of Computational Physics, 227, 3849-3877 [2] Restelli, M. & Giraldo, F.X. (2009). A conservative discontinuous Galerkin semi-implicit formulation for the Navier-Stokes equations in nonhydrostatic mesoscale modelling. SIAM J. Sci. Comp., 31, 2231-2257 [3] Kelly, J.F. & Giraldo, F.X. (2012). Continuous and discontinuous Galerkin methods for scalable nonhydrostatic atmospheric models: limited-area mode. Journal of Computational Physics, in review (2012). |
||||

12:30-13:30 | Lunch at Wolfson Court | |||

13:30-14:00 | Poster Session | |||

15:15-15:45 | Afternoon Tea | |||

17:00-17:25 | Evaggelos Kritsikis (Laboratoire des Sciences du Climat et l'Environnement (LSCE)) |
|||

Second-order conservative remapping between unstructured spherical meshes | Sem 1 | |||

Remapping from one finite-dimensional description of a function to another is a common problem in numerical modelling. In particular, information transit between meshes is required for, e.g., model coupling or mesh adaptation. In order to preserve the properties of a numerical scheme such as conservativity, accuracy, positivity, etc., the remapping algorithm must itself possess these properties. A second-order, conservative remapping between unstructured spherical meshes is presented. Areas are computed exactly by the defect formula and gradients estimated by the Gauss formula. Data is tree-structured, so that neighbour search is done in logarithmic time. In addition, the algorithm lends itself well to parallelisation. Numerical tests on various unstructured grids are given. |

Wednesday 26 September | ||||

Session: Vertical Session | ||||

09:00-09:25 | Yiyuan Li (Chinese Academy of Sciences) |
|||

The orthogonal curvilinear terrain-following coordinate for atmospheric models | Sem 1 | |||

This study designs a novel terrain-following coordinate, called the “Orthogonal curvilinear terrain-following coordinate” (Orthogonal sigma coordinate, OS-coordinate), which has the ability to handle both the well-known “high level errors” and the “pressure gradient force (PGF) errors” above the steep terrain in the classic terrain-following coordinate (sigma-coordinate). The OS-coordinate is proposed through a way quite different from and against the sigma-coordinate. The basis vectors of OS-coordinate which is unit and terrain-following are firstly designed, and then the corresponding definitions of every coordinates are solved, so that the scalar equations of OS-coordinate are solved. The PGF computational form has only one term in each momentum equation of OS-coordinate, thereby avoiding the PGF problem completely. Finally, a unified framework is proposed to combine the equations of the Cartesian coordinate, sigma-coordinate and OS-coordinate together via an “on-off”, and is used to implement idealized advection experiments for testing the OS-coordinate. The evaluation on the experimental results show the vertical coordinate surfaces in the OS-coordinate can be much smoother than those in the sigma-coordinate, and therefore the OS-coordinate outperforms all along the sigma-coordinate. This new coordinate can be used for oceanic models as well, and needs to be tested in both the idealized and real-case experiments. |
||||

09:25-09:50 | Takeshi Enomoto (Kyoto University) |
|||

Reexamination of non-hydrostatic formulations using the hydrostatic-pressure based co-ordinates | Sem 1 | |||

The rapid increase of computing power is making global non-hydrostatic simulations affordable. A natural approach is to extend the formulation to include the non-hydrostatic effect. The advantage of this approach is that the existing data assimilation and tools require minimal changes. ECMWF and JMA seem to pursue this approach. ECMWF has achieved TL7999 (corresponding to approximately 2.5 km) with a fast Lendre transform using the butterfly algorithm (Nils Wedi, pers. comm.). Hiromasa Yoshimura (MRI/JMA) has built a non-hydrostatic version of JMA GSM using double Fourier series. Their formulations are based on Laprise (1992) that proposes the vertical co-ordinates based on hydrostatic pressure. Juang (1992, 2000) also adopts hydrostatic sigma–co-ordinates in the vertical but there are subtle differences. The latter introduces the hydrostatic temperature. In a limited-area model, such as MSM, the hydrostatic temperature may be given externally. In a GCM, however, the hyd rostatic temperature must be determined internally if is not time-independent. We investigated the two formations and found the assumption of the hydrostatic state of Laprise (1992) may be used to diagnose the hydrostatic temperature within MSM. Similarly the hydrostatic assumption of Arakawa and Konor (2009) can be used. MSM is found to run stably with any of these diagnosed hydrostatic states. The diagnosed hydrostatic temperature would enable the application of the formulation of MSM to the global domain. |
||||

09:50-10:15 | Margarida Moragues (Barcelona Supercomputing Center (BSC-CNS)) |
|||

A Variational Multiscale Stabilized Finite Element Method to solve the Euler Equations for Nonhydrostatic Stratified Benchmarks | Sem 1 | |||

In this talk we present a Variational Multiscale Stabilization (VMS) for Compressible Euler Equations applied to the Finite Element (FE) solution of nonhydrostatic stratified flows. The VMS method was firstly presented by Hughes and co-workers [1] in the context of incompressible flows. In the present work, recentely presented in [3], we extend these concepts to Compressible Flows. In the framework of nonhydrostatic atmospheric dynamics, we test the algorithm for problems at low Mach numbers. A general version of the current compressible VMS technique was originally devised for Computational Fluid Dynamics (CFD) of compressible flows without stratification [2]. The present work is justified by the previously observed good performance of VMS and by the advantages that an element-based Galerkin formulation offers on massively parallel architectures, a challange for both CFD and Numerical Weather Prediction (NWP). Unphysical vertical oscillations that may appear for not well-balanced approximations are a relevant problem in NWP, especially in the proximity of steep topography. In that respect, to properly discretize the dominant hydrostatics, a particular interpolation technique is proposed. To evaluate the performance of the method in this context, some standard test cases of stratified environments are presented. References [1] T. Hughes, Multiscale Phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, CMAME 127 (1995) 387–401. [2] M. Moragues, M. V´azquez, G. Houzeaux, R. Aubry, Variational Multiscale Stabilization of Compressible flows in Parallel Architectures, Parallel CFD 2010, Taiwan, May 2010. [3] S. Marras, M. Moragues, M. V´azquez, O. Jorba, G. Houzeaux, A Variational Multiscale Stabilized Finite Element Method for the Solution of the Euler Equations of Nonhydrostatic Stratified Flows, J. Comput. Phys. (submitted 2012) |
||||

10:15-10:40 | Hiroe Yamazaki (University of Cambridge) |
|||

2D and 3D simulations of a nonhydrostatic atmospheric model on a block-structured Cartesian mesh | Sem 1 | |||

Under the rapid development of computing power, this study aims at developing a next-generation atmospheric model for ultra-high resolution simulations at horizontal grid intervals of less than 100 meters. Recently, Cartesian grids are drawing attention as an attractive choice for high-resolution atmospheric models that need to handle steep slopes in mountainous areas. They have the advantage of avoiding errors because of the slantwise orientation of grid lines in models using conventional terrain-following grids. In our model, a cut cell method is used for representing topography on a Cartesian grid. Also, a block-structured mesh approach is introduced to achieve computationally efficient Cartesian grid simulations with both high vertical resolution near the ground and reasonable conservation characteristics. The result of a 2D numerical simulation shows the model successfully reproduces a flow over a semi-circular mountain on a locally refined mesh around the mountain. The result agrees well with that using a uniformly fine mesh. Some recent 3D results of our model will also be discussed. |
||||

10:40-11:10 | Morning Coffee | |||

Session: Time Session | ||||

11:10-11:35 | Colm Clancy (Environment Canada) |
|||

Semi-implicit predictor-corrector methods for atmospheric models | Sem 1 | |||

A class of semi-implicit time integration schemes is proposed for use in atmospheric modelling. Various explicit predictor-corrector methods can be combined with an implicit treatment of the linear terms responsible for fast modes in order to improve computational stability. Linear analysis is used to identify promising algorithms. These are then tested in a finite volume shallow water model on an icosahedral grid. This framework has been developed at Environment Canada and will be extended to a more complete, three-dimensional atmospheric model in the future. Experiments with standard test cases show that the proposed time schemes allow for stable integrations with relatively long time-steps, while maintaining a sufficient level of accuracy. In particular, they are more robust than the traditional, single-stage semi-implicit method using the leapfrog discretisation and do not require a time filter to control the computational mode. |
||||

11:35-12:00 | Peter Lynch (University College Dublin) |
|||

Laplace transform integration and the slow equations | Sem 1 | |||

We consider the Laplace transform (LT) filtering integration scheme applied to the shallow water equations, and demonstrate how it can be formulated as a finite difference scheme in the time domain by analytical inversion of the transform. Both Eulerian and semi-Lagrangian versions of the scheme are analyzed. We show the relationship between the LT scheme and the slow equations. We demonstrate the advantages of the LT scheme by means of numerical integrations. |
||||

12:00-12:25 | Kevin Viner (Naval Research Laboratory) |
|||

A stable treatment of conservative thermodynamic variables for semi-implicit semi-Lagrangian dynamical cores | Sem 1 | |||

Atmospheric motions span a wide array of frequencies, the slowest of which provide us with our day to day weather. In numerical weather prediction it is necessary to narrow the focus to these lower frequencies for efficiency. The way this is typically done in global modeling is to apply an implicit time-differencing method to terms linked to high frequency motions while applying an explicit method to terms linked to low frequency motions like advection, thereby allowing a larger stable time step. While semi-Lagrangian schemes further increase efficiency by removing the stability restrictions of the standard CFL condition concerning velocity, they are still held to a slightly different deformation CFL condition concerning the local variation in velocity. As such, it is still necessary to slow these fast wave modes in the semi-Lagrangian framework to maintain a stable system. Semi-Lagrangian systems employing a conservative thermodynamic variable such as potential temperature as a prognostic variable face a unique dilemma in applying this standard method because the term responsible for gravity wave generation is also a vertical advection term which is absorbed into the total derivative. Application of the scheme in the absence of an explicit gravity wave term results in an unstable system since the fast gravity mode frequencies are not properly reduced. Stability can be maintained at the expense of both accuracy and efficiency by way of artificial damping and reduced time steps. The modification discussed here defines a basic state potential temperature field which is advected in an Eulerian fashion while the semi-Lagrangian method is applied to the perturbation potential temperature. As the gravity-wave term in question is now expressed explicitly through the total derivative of the basic state potential temperature, gravity mode stability is returned to the system; the semi-implicit treatment of the new term manages stability with respect to the associated CFL condition. Tests of the new scheme in the Navy Global Environmental Model (NAVGEM) show that the method allows stable integration at typical semi-Lagrangian time steps in the absence of artificial damping. |
||||

12:25-12:50 | Matthew Norman (Oak Ridge National Laboratory) |
|||

Multi-Moment ADER-Taylor Methods for Systems of Conservation Laws with Source Terms in One Dimension | Sem 1 | |||

A new integration method combining the ADER time discretization with a multi-moment finite-volume framework is introduced. ADER runtime is reduced by performing only one Cauchy-Kowalewski (C-K) procedure per cell per time step, and by using the Differential Transform Method for high-order derivatives. Three methods are implemented: (1) single-moment WENO [WENO], (2) two-moment Hermite WENO [HWENO], and (3) entirely local multi-moment [MM-Loc]. MM-Loc evolves all moments, sharing the locality of Galerkin methods yet with a constant time step during p -refinement. Five experiments validate the methods: (1) linear advection, (2) Burger's equation shock, (3) transient shallow-water (SW) , (4) steady-state SW simulation, and (5) SW shock. WENO and HWENO methods showed expected polynomial h -refinement convergence and successfully limited oscillations for shock experiments. MM-Loc showed expected polynomial h -refinement and exponential p -refinement convergence for linear advection and showed sub-exponential (yet super-polynomial) convergence with p -refinement in the SW case. HWENO accuracy was generally equal to or better than a five-moment MM-Loc scheme. MM-Loc was less accurate than RKDG at lower refinements, but with greater h - and p -convergence, RKDG accuracy is eventually surpassed. The ADER time integrator of MM-Loc also proved more accurate with p -refinement at a CFL of unity than a semi-discrete RK analog of MM-Loc. Being faster in serial and requiring less frequent inter-node communication than Galerkin methods, the ADER-based MM-Loc and HWENO schemes can be spatially refined and have the same runtime, making them a competitive option for further investigation. |
||||

12:50-13:30 | Lunch at Wolfson Court | |||

19:30-22:00 | Conference Dinner at St Catharine's College |

Thursday 27 September | ||||

Session: Tests Cases | ||||

09:00-09:25 | Todd Ringler (Los Alamos National Laboratory) |
|||

A Multi-Resolution Modeling Approach for Global Ocean Modeling | Sem 1 | |||

09:25-09:50 | Jared Whitehead (Los Alamos National Laboratory) |
|||

Potential Vorticity: A Diagnostic Tool for General Circulation Models | Sem 1 | |||

Maintaining correlation between tracers and the dynamical wind and temperature fields with which they interact is a desirable trait of climate and weather models. A systematic, explicit test is developed to measure the consistency between a dynamical core's integration of the momentum equation, and its tracer transport algorithm. Potential vorticity is used as a diagnostic tool allowing direct comparison between the treatment of the dynamics and the integration of passive tracers. Several quantitative and qualitative metrics are suggested to measure this consistency including grid-independent probability density functions. Comparisons between the four primary dynamical cores of the National Center for Atmospheric Research's (NCAR) Community Earth System Model 1.0 (CESM) are presented. It is found that the finite volume (CAM-FV) and spectral-element (CAM-SE) dynamical cores perform better than the spectral-transform based Eulerian (CAM-EUL) and semi-Lagrangian (CA M-SLD) dynamical cores in the presence of a breaking wave. |
||||

09:50-10:15 | Peter Lauritzen (National Center for Atmospheric Research) |
|||

A standard test case suite for two-dimensional linear transport on the sphere | Sem 1 | |||

It is the purpose of this paper to propose a standard test case suite for two-dimensional transport schemes on the sphere intended to be used for model development and facilitating scheme intercomparison. The test cases are designed to assess important aspects of accuracy in geophysical fluid dynamics such as numerical order of convergence, “minimal” resolution, the ability of the transport scheme to preserve filaments, transport “rough” distributions, and to preserve pre-existing functional relations between species/tracers under challenging flow conditions. The experiments are designed to be easy to set up. They are specified in terms of two analytical wind fields (one non-divergent and one divergent) and four analytical initial con- ditions (varying from smooth to discontinuous). Both conventional error norms as well as novel mixing and filament preservation diagnostics are used that are easy to implement. The experiments pose different challenges for the range of transport approaches from Lagrangian to Eulerian. The mixing and filament preservation diagnostics do not require an analytical/reference solution, which is in contrast to standard error norms where a “true” solution is needed. Results using a dozen state-of-the-art transport schemes that participated in the 2011 NCAR workshop on transport are presented. |
||||

10:15-10:40 | Mike Cullen (Met Office) |
|||

Evaluating numerical methods by using asymptotic limit solutions | Sem 1 | |||

Numerical models of the atmosphere and ocean 'solve' the governing Navier-Stokes equations, but because of the complexity of the true solutions, can only do so in a grossly averaged sense. Standard numerical analysis theory is then of limited use, because it expects solutions to be close to the truth. The Navier-Stokes equations are known to have solutions close to computable asymptotic limits in many cases of practical interest. This talk gives examples, including dynamics/boundary-layer interaction, frontal solutions where the asymptotic limit is singular, and long-term behaviour of baroclinic systems, together with numerical demonstrations. The technique should help to validate the new integration schemes being proposed. |
||||

10:40-11:10 | Morning Coffee | |||

Session: From Dynamical Cores to Operational Models | ||||

11:10-11:35 | Paul Ullrich (University of California, Davis) |
|||

MCore: Towards a Multi-Resolution Non-Hydrostatic Finite-Volume Dynamical Core | Sem 1 | |||

This talk will focus on recent efforts towards developing a true multi-resolution framework for modeling regional climate in a global domain. In particular, I will present ongoing work on the high-order non-hydrostatic finite-volume MCore atmospheric general circulation model, which has been designed for massively parallel systems and with a focus on using mesh refinement to seamlessly tie together regional and global scales. Further, efforts to develop test cases of intermediate complexity, which include moisture and simplified physics parameterizations, have been undertaken and simulated within MCore. These efforts will be presented, along with a look at some upcoming approaches to numerical modeling using finite-volume methods on the sphere. |
||||

11:35-12:00 | Jin Lee (National Oceanic and Atmosphere Administration) |
|||

A Three-Dimensional Finite-Volume Nonhydrostatic Icosahedral Modle (NIM) | Sem 1 | |||

The Nonhydrostatic Icosahedral Model (NIM) formulates the latest numerical innovation of the three-dimensional finite-volume control volume on the quasi-uniform icosahedral grid suitable for ultra-high resolution simulations. NIM’s modeling goal is to improve numerical accuracy for weather and climate simulations as well as to utilize the state-of-art computing architecture such as massive parallel CPUs and GPUs to deliver routine high-resolution forecasts in timely manner. NIM uses innovations in model formulation similar to its hydrostatic version of the Flow-following Icosahedral Model (FIM) developed by Earth System Research Laboratory (ESRL) which has been tested and accepted for future use by the National Weather Service as part of their operational global prediction ensemble. Innovations from the FIM used in the NIM include: * A local coordinate system remapped spherical surface to plane for numerical accuracy (Lee and MacDonald, 2009), * Grid points in a table-driven horizontal loop that allow any horizontal point sequence (A.E. MacDonald, et al., 2010), * Flux-Corrected Transport formulated on finite-volume operators to maintain conservative positive definite transport (J.-L, Lee, ET. Al., 2010), * All differentials evaluated as finite-volume integrals around the cells, *Icosahedral grid optimization (Wang and Lee, 2011) NIM extends the two-dimensional finite-volume operators used in FIM into the three-dimensional finite-volume solvers designed to improve pressure gradient calculation and orographic precipitation over complex terrain. The NIM dynamical core has been successfully verified with various non-hydrostatic benchmark test cases such as warm bubble, density current, internal gravity wave, and mountain waves. Physical parameterizations have been incorporated into the NIM dynamic core and successfully tested with multimonth aqua-planet simulations. |
||||

12:00-12:25 | Xingliang Li (China Meteorological Administration) |
|||

A multi-moment constrained finite volume model for non-hydrostatic atmospheric dynamics | Sem 1 | |||

Two dimensional non‐hydrostatic compressible dynamic cores for atmosphere are developed by using a new nodal type high order conservative method, so called multi‐moment constrained finite volume (MCV) method. Different from conventional finite volume method, the predicted variables (unknowns) in an MCV scheme are the values at the solution points distributed within each mesh cell. The time evolution equations to update the unknown point values are derived from a set of constraint conditions based on multi‐moment concept, where the constraint on the volume integrated average (VIA) for each mesh cell is cast into a flux form and thus guarantees rigorously the numerical conservation. Two important features makes MCV method particularly attractive as an accurate and practical numerical framework for atmospheric and oceanic modelling. 1) Using nodal values at uniformly located solution points as the predicted variables provides great convenience in de aling with complex geometry and source terms, and 2) High order schemes can be built by using constraints in terms of different moments, which makes the numerical schemes more flexible and efficient. We present in this paper the dynamic cores that use the third and the fourth order MCV schemes. We have verified the numerical outputs of both schemes by widely used standard benchmark tests and obtained competitive results. The present numerical cores provide a promising and practical framework for further development of non‐hydrostatic compressible atmospheric models. |
||||

13:35-14:00 | James Kent (University of Michigan) |
|||

DCMIP 2012: Tracer Transport Tests in Dynamical Cores | Sem 1 | |||

We investigate the accuracy of tracer transport algorithms in dynamical cores of weather and climate models using prescribed velocity test cases. The tests make use of three-dimensional flow, with the focus on extreme deformation, horizontal-vertical coupling, and flow in the presence of orography. These tests highlight common problems associated with advection in dynamical cores. Standard error measures allow us to assess convergence rates, monotonicity, and the effects of the diffusion mechanisms in tracer transport algorithms. We will present results form the 2012 Dynamical Core Model Intercomparison Project (DCMIP) workshop. |
||||

14:00-14:25 | Sasa Gabersek (Naval Research Laboratory) |
|||

Numerical Simulations with a Three-Dimensional Spectral Element Model | Sem 1 | |||

Highly accurate numerical methods for solving partial differential equations that have been traditionally used in the computational fluid dynamics have yet to be fully exploited for geophysical fluid dynamics applications, such as weather prediction. We will present results obtained with a fully compressible, non-hydrostatic spectral element (SE) model in three dimensions. All of our results are obtained using an eighth order polynomial (p=8), which provides the best compromise between accuracy and computational cost. The number of elements (h), into which the computational domain is decomposed, is varied among experiments to achieve different effective spatial resolutions. Introducing physical parameterizations into an SE model is challenging due to a number of aspects including the varying nodal spacing within each element. We highlight some of these challenges and our approaches to address them. The model is applied in a series of idealized experiments that include physical process parameterizations such as: i) flow over complex terrain with sub-grid scale mixing, and ii) a sea breeze that includes an evolving planetary boundary layer with surface fluxes. We discuss the issues that arise when physical parameterizations are introduced into SE models. Also, scalability results over many computational cores will be addressed. In addition, we will compare the efficiency and accuracy of the results with a finite difference model (FD). |
||||

14:25-14:50 | Celal Konor (Colorado State University) |
|||

Design of a dynamical core based on the nonhydrostatic unified system of equations | Sem 1 | |||

In this talk, we present the design of a dry dynamical core based on the nonhydrostatic “unified system” of equations. The unified system filters vertically propagating acoustic waves. The dynamical core predicts the potential temperature and horizontal momentum. It uses the predicted potential temperature to determine the quasi-hydrostatic components of the Exner pressure and density. The continuity equation is diagnostic (and used to determine vertical mass flux) because the time derivative of the quasi-hydrostatic density is obtained from the predicted potential temperature. The nonhydrostatic component of the Exner pressure is obtained from an elliptic equation. The main focus of this paper is on the integration procedure of this unique dynamical core. Height is used as the vertical coordinate, and the equations are vertically discretized on a Lorenz-type grid. Cartesian horizontal coordinates are used along with an Arakawa C-grid in the planar version of the dynamical core. The global version of the dynamical core is based on the vorticity and divergence predicting Z-grid formulation. The performance of the model in simulating a wide range of dynamical scales in the planar and global domains is demonstrated through idealized extratropical cyclogenesis simulations, and warm and cold bubble test cases. |
||||

14:50-15:15 | Almut Gassmann (Leibniz-Institut für Atmosphärenphysik, Kühlungsborn) |
|||

ICON-IAP: A non-hydrostatic global model designed for energetic consistency | Sem 1 | |||

The talk describes a new global non-hydrostatic dynamical core (ICON-IAP: Icosahedral Nonhydrostatic model at the Institute for Atmospheric Physics) on a hexagonal C-grid which is designed to conserve mass and energy. Energy conservation is achieved by discretising the antisymmetric Poisson bracket which mimics correct energy conversions between the different kinds of energy (kinetic, potential, internal). Because of the bracket structure this is even possible in a complicated numerical environment with (i) the occurrence of terrain-following coordinates with all the metric terms in it, (ii) the horizontal C-grid staggering on the Voronoi-mesh and the complications induced by the need for an acceptable stationary geostrophic mode, and (iii) the necessity for avoiding the Hollingsworth-instability. The model is equipped with a Smagorinsky-type nonlinear horizontal diffusion. The associated dissipative heating is accounted for by the application of the discrete product rule for derivatives. The time integration scheme is explicit in the horizontal and implicit in the vertical. In order to ensure energy conservation, the Exner pressure has to be off-centered in the vertical velocity equation and extrapolated in the horizontal velocity equation. Test simulations are performed for small scale and global scale flows. A test simulation of linear nonhydrostatic flow over a rough mountain range shows the theoretically expected gravity wave propagation. The baroclinic wave test is extended to 40 days in order to check the Lorenz energy cycle. The model exhibits excellent energy conservation properties even in this strongly nonlinear and dissipative case. The Held-Suarez test confirms the reliability of the model over even longer timescales. |
||||

15:45-16:10 | Nils Wedi (European Centre for Medium-Range Weather Forecasts (ECMWF)) |
|||

Global, non-hydrostatic, cloud-permitting, medium-range forecasts using the spectral transform method: progress and challenges | Sem 1 | |||

Numerical weather prediction (NWP) requires an answer in real time with a window of approximately one hour to run a medium-range global forecast such that it can be delivered in time to its customers worldwide. While computational efficiency remains one of the most pressing needs of NWP, it is an open question how to most efficiently use the computer power available over the next twenty years, while seeking the most accurate and cost-effective solution. At the same time, there are significant scientific challenges to increase resolution further, changes to the governing equations, and how sub-gridscale (SGS) processes are represented. ECMWF plans to implement a global horizontal resolution of approximately 10km by 2015 for its assimilation and deterministic forecast system, and approximately 20km for the ensemble prediction system (EPS). The scales resolved at these resolutions are still hydrostatic and the efficiency of the contemporary hydrostatic, semi-Lagrangian, semi-imp licit solution procedure using the spectral transform method is likely to remain a relevant benchmark. However, due to the relative cost increase of the Legendre transforms compared to the gridpoint computations, very high resolution spectral models may become prohibitively expensive. Moreover, spectral-to-gridpoint transformations require data-rich global communications at every timestep that may become too expensive on massively parallel computers. Recent progress in the development of fast spherical harmonics transforms (Tygert, 2008,2010) based on the butterfly scheme (O’Neil et al, 2010) mitigate the computational expense of the spectral transforms. Results are presented that demonstrate the cost-effectiveness of the fast Legendre transforms (FLT) on the ibm_power7 architecture. The FLTs save both memory and computing time enabling the "world's first" successful T7999 (or equivalently ~2.5km horizontal resolution) global weather forecast with a spectral transform model. Tygert, M., Fast algorithms for spherical harmonic expansions, II, J. of Comput. Physics, Vol. 227 (8), 2008, 4260-4279. Tygert, M., Fast algorithms for spherical harmonic expansions, III, J. of Comput. Physics, Vol. 229 (18), 2010, 6181-6192. O’Neil, M., F. Woolfe, V. Rohklin, An algorithm for the rapid evaluation of special function transforms, Appl. Comput. Harmon. Anal., Vol. 28(2), 2010, 203-226. |
||||

16:10-16:35 | Michail Diamantakis (European Centre for Medium-Range Weather Forecasts (ECMWF)) |
|||

Evaluation of mass fixing schemes for the ECMWF Integrated Forecasting System (IFS) | Sem 1 | |||

The purpose of this talk is to summarize the conservation properties of the ECMWF spectral semi-Lagrangian, semi-implicit IFS model and present results from some recently implemented algorithms for correcting global mass imbalances due to advection. Several established algorithms of different complexity, such as the ones described in [1,2,3,4], will be summarized and assessed in the context of NWP forecasts. Results from water species advection and passive tracer advection will be presented. Strengths and limitations of these schemes will be analyzed as well as their ability to preserve monotonicity and the impact they have on forecast skill. References Bermejo R., Conde J., 2002: A conservative Quasi-Monotone Semi-Lagrangian Scheme. MWR, 130, 423-430. A. Priestley, 1993: A Quasi-Conservative Version of the Semi-Lagrangian Advection Scheme, MWR, 121, 621-629. Zerroukat M., 2010, A simple mass conserving semi-Lagrangian scheme for transport problems, JCP, 229, 9011-9019. Williamson D.L., Rasch P.J., 1994, Water vapor transport in the NCAR CCM2, Tellus, 46A, 34-51. |
||||

16:35-17:00 | Michael Baldauf (Deutscher Wetterdienst (DWD)) |
|||

An exact analytical solution for gravity wave expansion of the compressible, non-hydrostatic Euler equations on the sphere | Sem 1 | |||

For the development and assessment of dynamical cores for atmospheric simulation models, suitable idealized test setups with known solutions are very useful. But only in rare cases exact analytical solutions exist for the underlying equation systems. In this work a slightly modified version of the original idea of Skamarock, Klemp (1994) is proposed: the quasi linear expansion of sound and gravity waves on a sphere induced by a weak warm bubble. For this case an exact analytical solution for the compressible, non-hydrostatic Euler equations was found for a shallow atmosphere and optionally with inclusion of Coriolis effects for a 'spherical f-plane-approximation'. This solution can be used as reference for convergence studies of global models. 'Small earth' convergence tests with the ICON model of the Deutscher Wetterdienst (DWD) and the Max-Planck Institut of Meteorology (MPI) are shown. |
||||

17:00-17:25 | Christiane Jablonowski (University of Michigan) |
|||

Highlights of the Dynamical Core Model Intercomparison Project (DCMIP) | Sem 1 | |||

The talk presents highlights of the Dynamical Core Model Intercomparison Project (DCMIP) and its associated 2-week summer school that has been held at NCAR in August 2012. DCMIP paid special attention to non-hydrostatic global models that are an emerging trend, especially in the climate modeling community. These allow high-resolution simulations and provide a pathway for embedded variable-resolution meshes in regions of interest. The primary goals of DCMIP were to assess the scientific designs of over 16 current and forthcoming dynamical cores, to establish new non-hydrostatic and moist dynamical core test cases in the dynamical core community, and to introduce innovative cyberinfrastructure tools. The talk surveys the DCMIP dynamical core test suite. This test hierarchy with increasing complexity (1) assesses the accuracy of 3D advection schemes via 3D deformational flows, (2) evaluates the evolution of gravity waves in non-rotating model configurations, (3) sheds light on the impact of orography, especially in the presence of orography-following vertical coordinates, (4) assesses dry baroclinic waves with dynamic tracers, and (5) suggests test cases of intermediate complexity that incorporate moisture and very simplified physical parameterizations. The latter include a new moist variant of the baroclinic wave test case, and idealized tropical cyclones. The test suite makes extensive use of small-planet experiments that are suggested as a computationally-efficient tool for high-resolution non-hydrostatic model evaluations. Highlights of the DCMIP model results are shown. |
||||

12:30-13:30 | Lunch at Wolfson Court | |||

15:15-15:45 | Afternoon Tea |

Friday 28 September | ||||

Session: From Dynamical Cores to Operational Models | ||||

09:00-09:25 | Günther Zängl (Deutscher Wetterdienst (DWD)) |
|||

The Icosahedral Nonhydrostatic (ICON) model: formulation of the dynamical core and physics-dynamics coupling | Sem 1 | |||

The nonhydrostatic dynamical core of the ICON model is formulated on an icosahedral-triangular C-grid and uses the edge-normal wind component, the vertical wind speed, density and virtual potential temperature as prognostic variables. In the vertical, a generalized terrain-following coordinate based on height is used that allows for a rapid decay with height of topographic structures, thereby reducing the numerical discretization errors over steep mountains. Moreover, a truly horizontal formulation of the horizontal pressure gradient term is used that greatly improves numerical stability over steep mountains. The spatial discretizations are second-order, with some slight degradation near the pentagon points of the basic icosahedron, and are primarily optimized for computational efficiency, ensuring strict mass conservation and consistent tracer transport but no exact conservation of energy or vorticity-related quantities. Time integration is based on a second-order predictor-corrector scheme that is fully explicit in the horizontal and implicit for the terms entering into vertical sound wave propagation. The dynamics time step is therefore limited by the CFL stability condition for horizontal sound wave propagation, but a longer time step (typically 4x or 5x) is used for tracer transport and fast-physics parameterizations. A important aspect of physics-dynamics coupling is that all terms related to latent heat release have to be converted into temperature changes at constant volume because density is used as a prognostic variable and is kept fixed during the physics call. |
||||

09:25-09:50 | Mark Taylor (Sandia National Laboratories) |
|||

CAM high-resolution AMIP simulations using a spectral finite element dynamical core | Sem 1 | |||

We will present results from high-resolution AMIP simulations using CAM-SE: The Community Atmosphere Model (CAM), running with the Spectral finite Element dynamical core from NCAR's High-Order Method Modeling Environment. CAM-SE is now supported within the Community Earth System Model. It uses fully unstructured conforming quadrilateral grids, and thus supports global high resolution through quasi-uniform grids and localized regions of high resolution with unstructured grids. For global 1/4 and 1/8 degree resolutions, CAM-SE runs efficiently on hundreds of thousands of processors on modern Cray and IBM supercomputers and obtains excellent simulation throughput. For variable resolution grids, our simulations benefit from a new mixed-formulation tensor-based hyper-viscosity operator. The mixed-formulation allows for efficient implementation in finite element methods, and a tensor form is used to provide the correct scale-aware dissipation, especially in the distorted elements which appear in resolution transition regions. |
||||

09:50-10:15 | John McGregor (CSIRO Marine and Atmospheric Research (CMAR)) |
|||

Comparing the formulations of CCAM and VCAM and their performance as atmospheric GCMs | Sem 1 | |||

Two cube-based atmospheric GCMs have been developed at CSIRO, the Conformal-Cubic Atmospheric Model (CCAM) and, more recently, the Variable Cubic Atmospheric Model (VCAM). The formulations of the dynamical cores of both models will be described and compared. CCAM is formulated on the conformal-cubic grid, whereas VCAM is cast on the equiangular gnomonic-cubic grid. CCAM is a 2-time-level semi-Lagrangian semi-implicit Eulerian model, whereas VCAM employs a split-explicit flux-conserving approach. Both models use reversible staggering for the wind components (McGregor, MWR, 2005) to produce good wave dispersion behavior. CCAM employs several orographic treatments that are not available for use in the VCAM dynamical core, with the interesting consequence that only VCAM requires hybrid vertical coordinates. Both models include the same efficient message-passing framework. Although VCAM avoids the message-passing overheads necessitated by the Helmholtz solver of CCAM, if does have some overheads from more frequent calls to the wind staggering/unstaggering routines. Aspects of the climatologies of both models will be compared and their overall advantages and disadvantages discussed. VCAM is presently being coupled to the Parallel Cubic Ocean Model (PCOM) from JAMSTEC; both the atmosphere and ocean model employ identical grids, and progress on this activity will be briefly presented. |
||||

10:15-10:40 | Thomas Melvin; Nigel Wood (Met Office) |
|||

The dynamical core of the Met Office's Unified Model | Sem 1 | |||

The Met Office's Unified Model (MetUM) is used across a wide range of spatial and temporal scales, from very high resolution nowcasting to centennial climate predictions. Key to the accuracy and efficiency of its predictions is the dynamical core. The scalability of the dynamical core is becoming increasingly critical to the efficiency of the model as a whole. An update on the progress, benefits and plans for the replacement dynamical core of the Met Office's Unified Model (ENDGame) will be given together with an overview of the project to design the next generation dynamical core for the the Unified Model (GungHo). |
||||

11:10-11:35 | Abdessamad Qaddouri (Environment Canada) |
|||

The CMC 15-Km operational deterministic global forecast system with Yin-Yang grid | Sem 1 | |||

At Canadian meteorological center (CMC), we are currently implementing the future 15-km operational deterministic global forecast system. In the horizontal we use spherical coordinates on the overset Yin-Yang grid while in the vertical, we use a log-hydrostatic-pressure coordinate on the Charney-Phillips grid. The global forecast on the Yin-Yang grid is performed by considering a domain decomposition (a two-way coupling method) between two limited-area models (LAM), discretized on the two panels of Yin-Yang grid and using the same time step. Each panel of the Yin-Yang grid system is extended by a static halo region that plays the same role as the piloting region in limited-area modeling. Since the two sub-grids of the Yin-Yang grid do not match, the update of the variables in the pilot region is done by interpolation. Each limited-area model uses the same fully implicit semi-Lagrangian method as in the GEM operational model to solve its own dynamic core. Geophysical fields are produced separately for each LAM grid however, The two LAM models use the same parametrization schemes for the physics configuration. For data assimilation, we use the Ensemble-Var (En-Var) approach where the 4D ensemble covariances are used to produce 4D analysis without TL/AD models. In this presentation we will present some preliminary results of this work. |
||||

11:35-12:00 | William Skamarock (National Center for Atmospheric Research) |
|||

The global nonhydrostatic atmospheric model MPAS: Preliminary results from variable-resolution mesh tests | Sem 1 | |||

The Model for Prediction Across Scales (MPAS) is comprised of geophysical fluid flow solvers using horizontally-unstructured spherical centriodal Voronoi meshes and a C-grid staggering of the prognostic variables. We have constructed a global atmospheric model for these meshes that solves the fully compressible nonhydrostatic equations using a finite-volume formulation coupled with a split-explicit time integration technique to handle acoustic modes. The Voronoi meshes are unstructured grids that permit variable horizontal resolution, and we plan to make use of variable-resolution meshes in our weather and regional climate applications. Towards this end, we have begun testing the robustness of variable-resolution global mesh simulations using both idealized test cases (baroclinic waves) and full NWP forecasts. Preliminary results for hydrostatic-scale tests (dx > 10 km) indicate that flow features are appropriately resolved relative to the local resolution if the mesh t ransition zones provide a smooth transition from coarser to finer resolution regions and if model filters are appropriately scaled by the local mesh spacing. We will present these results and discuss their implications for variable-resolution mesh configurations and for the end applications. |
||||

10:40-11:10 | Morning Coffee | |||

Session: Novel Optimisation Techniques | ||||

12:00-12:25 | Florian Prill (Deutscher Wetterdienst (DWD)) |
|||

The Icosahedral Nonhydrostatic (ICON) model: Scalability on Massively Parallel Computer Architectures | Sem 1 | |||

Simulation in numerical weather prediction and climate forecasting has a fast-growing demand for memory capacity and processing speed. For the last decade, however, computer technology has shifted towards multi-core chip designs while at the same time on-chip clock rates have increased only moderately. The parallel implementation of the ICON model's nonhydrostatic dynamical core therefore follows a hybrid distributed/shared memory approach, based on the Message Passing Interface (MPI) and the OpenMP API. The ICON code couples the different encapsulated components of the earth system model, e.g. dynamics, soil, radiation, and ocean, with high-level language constructs. Its communication characteristics and programming patterns are designed to meet the main challenges in high performance computing, i.e. load balancing, cache efficiency, and low-latency networking, and take the unstructured triangular C-grid into account, which implies indirect addressing. Besides basic optimization strategies such as loop tiling and grid point reordering, the implementation employs special domain decomposition heuristics, parallel range-searching algorithms with logarithmic complexity, and makes use of asynchronous I/O servers to deal with the potentially prohibitive amount of data generated by earth system models. This facilitates the ICON code to extract an adequate level of performance on a wide range of HPC platforms, targeting large scalar cluster systems with thousands of cores as well as vector computers. |
||||

13:30-13:55 | Hirofumi Tomita (RIKEN/AICS) |
|||

SCALE-LES: Strategic development of large eddy simulation suitable to the future HPC | Sem 1 | |||

The Large Eddy Simulation is a vital dynamical framework to investigate the cloud-aerosol-chemistry-radiation interaction from the viewpoint of climate problem. So far, the LES using in the meteorological filed are having several problems. One problem was that it is large grid-size used, compromising to the suitability of LES. In addition, the aspect ratio of horizontal and vertical grids was much larger than unity. The grid-size must be reduced to several 10m and it is desirable that the aspect ratio is near unity for the atmospheric LES. The target domain was also narrow for less of computer resources. The large-scale computing using the recent powerful super-computer may enable us to conduct the LES with reasonable grid-size and wide domain. Ultimately, the global LES is one of milestones in near future. Another problem in LES applied on meteorological field is that the heat source owing to water condensation is injected in a grid box. Strictly considering, the grid-box heating collapse the theory of LES that the grid size is in the energy cascade domain. Nevertheless, we have used the dry theory of LES. Beside the above problem that should be resolved in the future, we are now confronting with computational problems for such large-scale calculations. The numerical method of fluid dynamical part in the atmospheric model has been shifted from the spectral transform method to the grid-point method. The former is no longer acceptable on the massively parallel platforms form the limitation of inner-connect communication. On the other hand, the latter also contains a new problem, which is so-called memory bandwidth problem. For example, even on K Computer, the B/F ratio is just 0.5. The key to get high computational performance is the reduction of load/store from and to the main memory and efficient use of cash memory. Similar problem occurs in the communication between computer nodes. The multidisciplinary team (Team SCALE) in RIKE/AICS is now tackling to such problems. |
||||

13:55-14:20 | Ryuji Yoshida (RIKEN/AICS) |
|||

An Inter-Comparison of Icosahedral Climate Models on the G8 Call: ICOMEX Project | Sem 1 | |||

The ICOsahedral-grid Models for EXascale Earth system simulations (ICOMEX) is the consortium for the climate models development toward the exascale computing. It started since October 2011 as one of the G8 call projects. Participated are NICAM (Japan), ICON (Germany), MPAS (UK and US), and DYNAMICO (France) model teams. On the road to the exascale computing, we would find many road blocks, such as file I/O speed, lower byte/flops rate, outer/inner node communication. In order to examine these problems, six working groups are seated together. The Japan team works on the model inter-comparison to be synergistic among working groups. Although all the participated models use the icosahedral-grid, the discretization methods are different. For example, the hexagonal shape of a control volume is used in NICAM, while the triangular shape is used in ICON. Through intercomparison of both computational and physical performances between models, we will find which aspects of the model configuration are advantagous toward exascale computing. Two types of experiments were performed until now using NICAM and ICON: the baroclinic wave test (Jablonowski and Williamson, 2006) and the statistical climatology test (Held and Suarez, 1994). The horizontal resolution is from 240km (glevel-5) to 14km (glevel-9), and higher resolution runs will be tested in future. For the Held and Suarez test case, NICAM and ICON simulated climatology similar to that shown in the original paper. The baroclinic wave test is also compared between the two models. To investigate computational aspects, we examined strong scaling of parallel computing. Tests were performed on the Westmere and Bulldozer machine by using 5 through 40 MPI processes. We found that the two models have a good scaling: the measured scaling efficiency is 0.8-0.9. We plan to perform the intercomparison experiments using K computer using a larger number of processes O(10^5) to examine detail profiles of very massive parallel cores. |
||||

14:20-14:45 | Jean Côté (UQAM - Université du Québec à Montréal) |
|||

Time-parallel algorithms for weather prediction and climate simulation | Sem 1 | |||

The forecast of weather relies on computer models that need to be executed in real-time, meaning that a forecast needs to be disseminated to users well before the time period for which it is made. A challenge in the future will be to succeed in using the computing power available in massively parallel high-performance computers and meet the real-time requirement. Until now weather forecast and related climate simulation models have taken advantage of the parallelism of the computers by dividing the task to be performed in the horizontal space dimensions. The purpose of this work is to develop algorithms that allow also parallelism in the time dimension. This increased parallelism should allow an acceleration of the execution time of weather and climate models. This acceleration in turn permits an increase in the space accuracy of models while still meeting the real-time requirement. The talk will present our preliminary work on the "Parareal" algorithm that has been developed for that purpose (Lions et al., 2001) and whose applications to date have included among others air quality, but ignored weather forecasting. Weather forecasting presents a challenge for the method because of the presence of waves and advection. An important question is to examine how the traditional way to accelerate models with the semi-implicit semi-Lagrangian methodology can be advantageously blended with the “Parareal” approach. |
||||

12:30-13:30 | Lunch at Wolfson Court | |||

14:45-15:15 | Afternoon Tea |