# THEMATIC PROGRAMS

November 30, 2015
 Takayuki Aoki Qiang Du David A. Randall Jorn Behrens Max Gunzburger Todd D. Ringler G. Beylkin Daniel Guo Thomas E. Rosmond Luca Bonaventura Thomas Heinze Masaki Satoh V. Cheruvu Detlef Hinneburg Masaki Satoh 2 Elisabetta Cordero Mariano Hortal R.K. Scott J. Coté 1 Christiane Jablonowski M. Shoucri J. Coté 2 Lili Ju Andrew Staniforth J. Coté 3 Eigil Kaas Reiji Suda Terry Davies Oswald Knoth S. Sumetkijakan John Drake S.J. Lin Yukio Tanaka A. Fournier Bennert Machenhauer Hirofumi Tomita M. Fox-Rabinovitz L.M. Polvani Motohiko Tsugawa Francis X. Giraldo 1 Janusz Pudykiewicz Agathe Untch Francis X. Giraldo 2 Abdessamad Qaddouri 1 K.S. Yeh, Koji Goto Abdessamad Qaddouri 2 Seong Young Yoon Hiroshi Yoshida

Takayuki Aoki, Tokyo Institute of Technology and Hiroshi Yoshida, Tokyo Institute of Technology
Shallow water equations in spherical geometry solved by a high-accurate IDO scheme and overset grid

A high-accurate IDO (Interpolated Differential Operator) Scheme[1] has been applied to shallow water equations in spherical geometry represented by Williamson's test cases[2]. The IDO scheme required not only the discretized physical value on the grid point, but also its spatial derivatives as independent variables. The additional variables require solving more equations derived from the given equations, so that it makes possible to construct a high-accurate interpolation function around the local area of the grid point. The interpolation function reduces to an approximation solution of the given partial differential equation. For Williamson's test case 2, our result keeps the initial balance much better than 4th-order finite difference scheme. The result of the test case 5 has good agreement with that of Pseudo-Spectral method.
In order to avoid the singularities at the poles, we introduce two overset grids to replace the area around the poles. For the communication between the grid boundaries, third-order interpolated is used due to making use of spatial derivatives. Rossby-Haurwitz wave of the test case 6 propagates across the grid overset region without any oscillations, and is in good agreement with the case of a single grid. Any problems will not expected to apply to massively parallel computers.

References

[1] T. Aoki, "Interpolated Differential Operator (IDO) Scheme for Solving
Partial Differential Equations", Comput. Phys. Comm. Vol.102, No.1-3
(1997) 132-146.
[2] Williamson et al., "A Standard Test Set fir Numerical Approximations
to the shallow", J. Comp. Phys., Vol.102 (1992) 211-224.

Jorn Behrens, Munich University of Technology
Achieving Mass Conservation in Adaptive Semi-Lagrangian Advection Schemes

Semi-Lagrangian methods receive wide acceptance for their favorable properties, namely stability, efficiency, parallelization. Even adaptivity with unstructured, locally refined meshes is easily achieved with this method of discretization.

However, semi-Lagrangian methods in their simple form lack conservation properties. In many fields of application, e.g. climate modelling, or modelling of nonlinear phenomena with shocks, conservation of mass (and other physical constituents) is essential.

We propose several modifications which add mass conservation to the advection case of the semi-Lagrangian method. A comparison highlights properties of the proposed schemes in different situations. Applications in planar and spherical geometries are given. The effectiveness of mass conservation is shown in test cases with diverging and converging wind fields.

Luca Bonaventura, Max Planck Institute for Meteorology Hamburg
Project for the development of a nonhydrostatic dynamical core for climate and NWP simulation.

The joint project of Max Planck Institut for Meteorology Hamburg and Deutscher Wetterdienst to develop a new dynamical core for climate and weather forecasting simulation is presented and outlined. The project has started in the spring of 2002 at MPI Hamburg and is expected to produce a new dynamical core by the end of 2005.

A dynamical kernel based on the discretization of the fully compressible, non-hydrostatic equations is under development. The numerical model will use a icosahedral grid for the horizontal discretization and employ a height based vertical coordinate. Cell centered finite volume discretizations of the continuity equation on either the hexagonal cells or the dual triangular cells are going to be employed.

As a first step towards the full model, conservative, semi-implicit numerical schemes to discretize the shallow water equations the sphere are being developed. Accuracy and benefits of various possible formulations will be analyzed by comparing results obtained on standard test cases.

Preliminary results of numerical simulations will be presented.

Elisabetta Cordero, Met Office
Normal mode analysis of the New Dynamics'' discretisation employed at the Met Office

The set of equations governing the atmospheric motions employed in the operational Met Office NWP and global climate model is the non-hydrostatic and deep-atmosphere form of the Navier-Stokes equation set. The equations are discretised on a uniform lat-long Arakawa C-grid in the horizontal and a Charney-Phillips grid in the vertical using a semi-implicit semi-Lagrangian finite difference time integration (New Dynamics'').

A key aspect in the development and testing of the New Dynamics was the creation of 1D and 2D versions of the model. These versions were used to investigate the numerical stability of the scheme, and to test its behaviour in the presence of orography and for different top and bottom boundary conditions. Some results from these idealised models have been presented at previous PDEs on the sphere''
workshops.

In order to bridge the gap between these experiments and various linear analyses, the linear normal modes of the exact 2D New Dynamics discretisation have been computed. The program has the flexibility to be used to analyse the 1D discretisation and different vertical staggering of the variables, and includes switches to change the top and bottom boundary conditions.

At present it has been used to compute the normal modes of the 1D column model with both an isentropic and an isothermal basic state. In the isentropic basic state case, the numerical results are found to agree well with the analytic ones. The normal modes of the isothermal basic state 1D column model have been used to initialise the 1D New Dynamics: some of the numerical results will be presented.

Future work will include: completing the analysis of the properties of the 1D and 2D versions of the existing discretisation, for different boundary conditions and vertical staggering of the variables; using the established normal mode model as a tool to analyse the properties of improved alternative discretisations, before implementing them in the full NWP model; and extending the program to include the trajectory calculation.

Terry Davies, Met Office, Bracknell, UK
Three dimensional test problems for global atmospheric models

There are a wide range of idealised tests available in 1 or 2 dimensions for testing atmospheric code but only limited 3 dimensional tests. The dynamical core tests of Held and Suarez need to be run for several hundred days to approach quasi-equilibrium. Full 3 dimensional balanced flows can be easily set up on an f-plane in limited area models where the lateral boundaries can be used to maintain a constant (in time) large scale flow. Here I will describe a similar approach for global tests where the flow is close to geostrophic and hydrostatic balance when orography is absent. Orography can be introduced inducing perturbations to the flow which increase with either the size of the orography or the speed of the flow. These tests are 3 dimensional analogues to some of the shallow-water test set and are useful in testing dynamics algorithms and the response to orographic forcing. The tests need to be run for no more than 10 days to provide meaningful results making them ideal for testing at high resolution.

John Drake and Daniel Guo
Smooth Grid Transformations with Spectral Methods for the Shallow Water Equations

The use of the Schmidt transformation in a semi-implicit, semi-Lagrangian formulation for shallow water equations in spherical geometry gives the ability to focus resolution over one particular area of interest. The spectral transform computation is modified, to account for the grid transformation, by solving a linear system in spectral space. This paper will present a formlation for advection of vorticity and divergence that fits well with smooth grid transformations. Numerical results using this formulation for the shallow water test cases will be presented and generalizations to other types of grid transformations will be discussed.

A. Fournier, University of Maryland, G. Beylkin, University of Colorado, V. Cheruvu, University of Colorado, S. Sumetkijakan, University of Maryland
Multiresolution Adaptive Spectral Elements: Application to shallow-water flow on the sphere

We present a high-order, adaptive numerical method for numerically solving advection-diffusion PDEs, and demonstrate the advantage of this method for fluid dynamics problems, especially that generate quasi-singular or coherent structures.

There are four aspects of the method that we will elucidate. First, the spatial discretization employs piecewise-polynomial spaces similar to
the traditional spectral-element method; however we use orthogonal multiresolution analysis (MRA) comprising many nonconforming elements of various scales. This MRA is constructed using multiwavelets, which are orthonormal basis functions whose span efficiently describes the extra detail'' information gained in refining some of the spectral elements (Alpert et al. 1999).

Second, the adaptivity criterion for a set of elements is the norm of the projection onto the corresponding multiwavelets. Using this simple
criterion, the algorithm refines and coarsens the representation of steepening or translating structures. Our criterion replaces more traditional criteria, based on local gradients or local spectra.

Third, we take advantage of the spatial discretization, to improve the time-stepping scheme. We use the exact linear part'' scheme (Beylkin et al. 1998), which amounts to using scaling and squaring to compute the exponential propagator due to the diffusive linear operator. This adaptive operator exponentiation is made feasible by our spatial discretization. It significantly improves convergence for high Reynolds numbers and large time steps.

Finally, we use block-sparse data structures to achieve efficiency.

We will demonstrate these features using the standard shallow-water test suite of Williamson et al. (1992).

Francis X. Giraldo, Naval Research Laboratory
A Spectral Element Semi-Lagrangian Atmospheric Model (SESLAM)

A new dynamical core for NWP based on the spectral element semi-Lagrangian method developed in [1] and [2] is presented.
This work extends the SEE-AM dynamical core (also presented at this workshop) to a semi-Lagrangian time discretization, and to a spectral element vertical discretization. The semi-Lagrangian time discretization will allow for much larger time-steps while the spectral element vertical discretization will cure many of the ailments of our current finite difference method. In addition, this new vertical discretization will allow the SESL scheme to be applied to all three spatial dimensions. A progress report on this effort will be given including results for at least two test cases: the Rossby-Haurwitz waves 1 and 4.

[1] F.X. Giraldo, The Lagrange-Galerkin spectral element method on
unstructured quadrilateral grids, Journal of Computational Physics, Vol 147,
114-146 (1998).

[2] F.X. Giraldo and J.B. Perot, A spectral element semi-Lagrangian (SESL)
method for the spherical shallow water equations,
Journal of Computational Physics, in review (October 2001).

Francis X. Giraldo, Naval Research Laboratory and Thomas E. Rosmond
A Scalable Spectral Element Eulerian Atmospheric Model (SEE-AM): Dynamical Core Tests

A new dynamical core for NWP based on the spectral element method is presented. In this work, the 3D primitive hydrostatic atmospheric equations are written, discretized, and solve in 3D Cartesian space. The advantages of this approach are: the pole singularity which plagues all gridpoint methods disappears, the horizontal operators can be approximated by local high-order elements, and any grid can be used including lat-lon, icosahedral, hexahedral, and adaptive unstructured grids. The locality property of spectral elements means that the method will scale efficiently on distributed-memory computers. In order to validate our 3D atmospheric dynamical core,
we have run three test cases: Rossby-Haurwitz waves 1 and 4, and the Held-Suarez test case. Comparisons with the Navy's operational NWP model (NOGAPS) using the Rossby-Haurwitz waves demonstrate the high-order accuracy of the solutions obtained with our new model. The model is shown to scale quite well on distributed-memory computers.

[1] F.X. Giraldo, A spectral element shallow water model on spherical geodesic grids, International Journal for Numerical Methods in Fluids, Vol. 35, 869-901 (2001).

Max Gunzburger, Iowa State University, Lili Ju, Iowa State University and Qiang Du, Penn State University
Solution of PDEs on the sphere by finite volume and finite element methods using CVT grids

Centroidal Voronoi tessellations (CVTs) are special Voronoi tessellations for which the generators are also the centers of mass (with respect to a given density function) of the the corresponding Voronoi regions. Such tessellations have many remarkable properties and are especially well suited for efficiently generating very high-quality grids. We discuss CVTs on the sphere and their properties, and also present some very efficient, eminently parallelizable, probabilistic algorithms for their construction. We demonstrate the superior uniformity of "uniform" CVT grids and also how locally refined grids can be systematically constructed within the CVT framework. We then study, both theoretically and computationally, the use of uniform and locally refined CVT grids for the solution of model PDE problems on the sphere by finite volume and finite element discretization methods. Among the result we present are optimal error estimates for the approximate solutions and computational results which demonstrate that optimally accurate results can be obtained through the use of both uniform and nonuniform CVT grids. For example, for finite volume discretizations, CVT grids yield second-order accuracy (with respect to L2 norms).

Thomas Heinze, Technische Universitaet at Muenchen, Germany
A close look on test case 5

In 1988 L.L.Takacz introduced the flow over an isolated mountain as a test problem for the shallow water equations (SWE). After Williamson et.al. added it as fifth test case to the NCAR test suite for the SWE many researchers investigated this problem.

One of the big challenges for the next generation of global circulation models (GCM) will be adaptive meshes. The most suitable test case for this kind of grid generation out of the test suite is test case 5. Because of this it will become more and more important to understand the inherent problems of test case 5.

We investigate test case 5 with three different models (GME, SEAM, FEMmE) and compare the results. Depending on the model we observe different problems. After the discussion of them we will propose an additional test and a modified mountain to get a deeper understanding of the model numerics.

Christiane Jablonowski, University of Michigan
New idealized test cases for dynamical cores

Idealized test cases for primitive equation based dynamical cores have become a standard assessment tool for the strengths and weaknesses of current dynamics packages in GCMs. But although tests like the Held-Suarez test or the breaking polar vortex pattern show important model characteristics, additional test cases will be needed for future model intercomparisons.

The talk will propose a new standard test sequence for dynamical cores with pressure-based vertical coordinates. The test suite comprises two parts. First, the model will be initialized with steady state, balanced initial conditions that are an analytical solution to the primitive equations. This initial flow consists of a zonally symmetric basic state with two jets in midlatitudes and a realistic temperature profile. The careful design guarantees static, inertial and symmetric stability properties, but is unstable with respect to baroclinic or barotropic instability mechanisms. Model integrations over thirty days reveal how well the model can keep its initial state. Second, a well-resolved small amplitude Gaussian hill perturbation is superimposed on the initial state. This triggers baroclinic wave activities that lead to explosive cyclogenesis.

The test sequence has been applied to NCAR's newest dynamical cores (CAM2 model framework), to NASA's finite volume model (developed at GSFC) and the icosahedral model GME (developed at DWD). In addition, model results based on cubed meshes and reduced grids are available. The intercomparison reveals interesting characteristics that have already led to model improvements. The diagnostics include error norm statistics, grid point data comparisons and wave number analyses.

In the future the test will serve as a standard test case for adaptive grid simulations that are capable of tracking localized phenomena at high resolutions. A brief outline of this project will be given.

Eigil Kaas, Danish Meteorological Institute and Bennert Machenhauer
An accurate cell-integrated semi Lagrangian and semi-implicit scheme based on step-functions

A new and very accurate cell-integrated semi Lagrangian (CISL) two time level scheme has been formulated and tested for the shallow water equations in a plane channel model with realistic variation of the Coriolis-parameter and including topography. The implementation includes semi-implicit time stepping of the gravity wave terms. Furthermore, two-dimensional advection of passive tracers has been tested in idealized flows as well as in a fully general flow.

Regarding basic features, the new formulation is based on step-functions, which implicitly are advected with the flow across cell boundaries. The scheme is exact in case of passive advection by a flow that is constant in time and space. This means that the accuracy is of indefinite order. In addition to this attractive feature the scheme shows very small numerical dispersion in fully non-linear and divergent flows and in such flows the scheme is locally (and globally) mass conserving, positive definite and monotonic. The high order of accuracy is achieved by introducing two additional variables. For a given prognostic variable the memory carrying variables are: the grid-cell integrals, the values at the grid cell corner points and the functional average on the intersections between grid cell corner points in either of the two spatial directions. The penalty of running the new scheme is an increase in memory consumption by a factor of three.

In the present preliminary formulation of the scheme only the mass field is treated with the CISL scheme, while a traditional semi-Lagrangian scheme is used for the momentum equations. It is, however, possible to use the CISL scheme also for the cell-integrated form of the momentum equations. In this case also momentum is conserved.

Results of test integrations will be presented and compared to integrations based on finite difference and spectral formulations of traditional semi-implicit Eulerian and semi Lagrangian shallow water models. The very convincing tests of passive advection by idealized flows and by a fully divergent flow simulated by the shallow water model are also shown.

The presentation is completed by a discussion of the generalization of the scheme to fully three-dimensional flows on the sphere. Briefly, the semi-implicit implementation follows Machenhauer and Olk (1997) and the application on the sphere will follow the formulation presented in Nair and Machenhauer (2002) with some modifications. Following Machenhauer and Olk (1998) the generalization to three dimensions is based on two-dimensional CISL advection on model levels in combination with a diagnosing of the vertical advection under the assumption of hydrostatic balance. This formulation ensures a more consistent treatment of the vertical advection problem than is the case in the traditional hydrostatic semi-Lagrangian models based on three dimensional trajectories.

REFERENCES
Machenhauer, B. and M. Olk, 1997: The implementation of the semi-implicit scheme in cell-integrated semi-Lagrangian models. The Andre J. Robert Memorial Volume. Numerical Methods in Atmospheric and Oceanic Modelling. (C. Lin, R. Laprise, H. Ritchie, Eds.) CMOS /NRC Research Press, Ottava, Canada. 103-126.

Machenhauer, B. and M. Olk, 1998: Design of a semi-implicit cell-integrated semi-Lagrangian model. MPI Workshop on Conservative
Transport Schemes, MPI for Meteorology Report No. 265. 76-85.

Nair, R. and B. Machenhauer, 2002: The Mass-Conservative Cell-Integrated Semi-Lagrangian Advection Scheme on the Sphere, Mon. Wea. Rev., Vol 130, No.3, 649-667.

Oswald Knoth, Institute for Tropospheric Research and Detlef Hinneburg, Institute for Tropospheric Research
A nonhydrostatic anelastic model with a cut cell approach and implicit time stepping

The implementation of a global nonhydrostatic anelastic model on the sphere in a lat-lon-z grid is presented. The representation of the orography is realized by cut cells, which describe approximately the intersection of the orography boundary with the grid cells. Therefore the free fluid part in each grid cell is characterized by the free cell volume and the free face area of the six cell faces. Since the lat-lon grid is locally orthogonal the discretization in space can be expressed with the above defined cell characteristics.

The spatial discretized system is solved adaptive in time by a combination of a Rosenbrock-method for the advcetion-diffusion part and a Chorin-type projection method for the pressure. This implicit time integration procedure avoids time step restrictions which are caused by small volume cells at the poles and at the cutting boundary and by fast but unimportant physical processes. The resulting linear systems are solved by preconditioned CG-like methods. For the advection-diffusion part the BICGStab method is used. The preconditioner is built up by an approximate matrix factorization of the transport terms (advection, diffusion) and the source terms (Coriolis, curvature, buoyancy).
The preconditioner for the pressure equation is a multigrid method with a plane smoother to overcome the anisotropies in the grid due to the large horizontal/vertical grid ratio and the pole singularity.

The code is parallelized by a decomposition of the computational domain in rectangular blocks in horizontal as in vertical direction. These blocks are distributed than on the available processors. The decomposition allows also to apply different grid resolution in the individual blocks. One possible application is the use of a coarser resolution in the polar regions. The parallelization requires changes in the solvers for the linear systems. Especially for the pressure computation we use an extended pressure equation with additional flux variables (pressure gradients) at the boundaries between neighbouring blocks. We will present computational results for simplified test scenarios.

L.M. Polvani, Princeton University and R.K. Scott, Columbia University
An initial-value test case for dynamical cores of atmospheric general circulation models

The development of general circulation models (GCMs) is an important task among the current efforts to understand and predict the climate. A key element of atmospheric GCMs is the so-called "dynamical core", the component of a GCM that deals with the numerical solution of the dry, adiabatic primitive equations. In order to build reliable dynamical cores, it is important to test them against known solutions.

At present, the only widely used test case for dry dynamical cores is the one proposed by Held & Suarez (1994), in which simply defined parametrizations of two key physical processes (thermal relaxation and surface drag) are added to the dynamical core. The model behavior is then tested by performing a 1,000-day integration and comparing the time-averaged fields to those in the Held & Suarez test case. One major drawback of this test case, however, is that a large amount of averaging needs to be performed before the results of a given model can be compared to those of the test case. It is conceivable that subtle coding errors or noisy features may not reveal themselves owing to the averaging. A second drawback is that numerical convergence was not demonstrated by the authors. A final, though minor drawback, is that new parameterizations need to be added to the dynamical core itself.

We here propose a new test case, meant to address the shortcomings of the Held & Suarez test case, and serve as a complementary tool for testing dynamical cores. Our new test case is an initial-value problem, with the zonal winds specified to be those of a baroclinically unstable, midlatitude zonal jet, analytically specified to be very close to the LC1 paradigm described by Thorncroft et al (1993). This zonal jet is slightly perturbed, and its evolution is integrated for 10 days. We present snapshots of the fields at various time intervals (e.g. the vorticity near the surface) as the baroclinic instability develops. We also show the time evolution of several diagnostic quantities over the 10 days of integration, and provide tables of these to serve as numerical benchmarks against which new dynamical cores can be precisely and quantitatively compared.

Unlike the Held & Suarez test case, our test case involves the addition of no additional parametrizations to the dynamical core, and requires a much shorter integration time. Also, the instantaneous fields computed with a new dynamical core can be compared directly to those of the benchmark, without the need for temporal or spatial averaging. Finally, by demonstrating numerical convergence, our new test case de facto provides a new, non-trivial, exact -- albeit numerically derived -- solution to the time-dependent primitive equations in spherical coordinates.

Janusz Pudykiewicz, Meteorological Service of Canada
Numerical simulation of reactive flow in spherical geometry

The set of equations governing reactive flow in spherical geometry describes a complex system in which the fluid dynamics is coupled with local processes such as chemical reactions or phase changes. In mathematical terms the reactive flow equations consist of Navier-Stokes equations coupled with continuity equations for a number of interacting scalar fields characterizing the chemical composition of fluid. In order to solve this very complex set of equations in an efficient manner we first discretise the spatial derivatives on an icosahedral mesh defined in Cartesian coordinates. The system of ordinary differential equations obtained following this procedure is then solved using several different numerical methods including the third and fourth order Runge-Kutta schemes. A discussion of constraints which are required in order to maintain the monotonicity of the the scheme will be also presented. In particular, it will be shown that the solver is both mass conserving and nonoscillatory which make it ideally suited for the solution of the complex reactive flow problems in spherical geometry. This property will be illustrated by the examples of application of the scheme for the simulation of tropospheric chemistry. In conclusion, the advantages of using geodesic grids for the simulation of reactive flows will be summarized.

Parallel Elliptic Solvers for the GEM Model

Numerical weather prediction using the Canadian GEM model involves the solution of separable elliptic boundary value (EBV) problems in spherical geometry. This type of equation arises either as Helmholtz or as horizontal diffusion problems.

The direct solution of the EBV problem involves a transform, in the variable-mesh case a full-matrix multiplication, where the cost per grid point rises linearly with the number of grid points along the transform direction. In order to improve the performance of the EBV problem, the matrix product in the direct solution is accelerated by using either the Strassen method or by exploiting the symmetry of the mesh. An iterative solution of the EBV has also been implemented and compared to the direct solution.

The purpose of this presentation is to report on the details of these implementations and to show the improved performance of the EBV problem solution either by accelerating the full-matrix in the direct solver or by using a preconditioned Conjugate Gradient iterative solver. The direct and iterative solvers have been tested on the NEC SX-4 and SX-5.

Todd D. Ringler, Colorado State University and David A. Randall, Colorado State University
Impact of an Energy Conserving Scheme in the CSU Geodesic-Grid AGCM

This talk will review a new dynamical formulation of the governing equations that has recently been implemented into the CSU AGCM. The new dynamical formulation is characterized by its conservation of mass, tracers, and total energy. Furthermore, the formulation exhibits a consistency between the momentum form of the governing equations and the vorticity-divergence form of the equations.

We will compare simulations that use this new dynamics package to simulations that use an older formulation. We will compare and contrast the results in terms of energy and enstrophy spectra. Furthermore, using a full physics AGCM, we will show that the new formulation leads to an changes not only in the dynamical fields, but also in the parameterized fields such a precipitation and cloud cover.

We have extended the new dynamical scheme to include the discretization of second order tensors by extending the definition of the discrete gradient operator to operate on vectors, as well as scalars. In addition, we have developed a conservative form of the divergenceof second order tensors. We will present results on the impact of this discretization in simulations of the 2-D turbulence.

Masaki Satoh, Frontier Research System for Global Change/ Saitama Institute of Technology, Japan
Conservative scheme for a non-hydrostatic climate model

For the equation set of the non-hydrostatic icosahedral global climate model being developed at Frontier Research System for Global Change, we have devised a conservative scheme of the compressible non-hydrostatic equations. The scheme is based on the flux form equations of total density, three components of momentum, total energy, and densities of water substances. Time-splitting is used with the leap-frog or 2nd order Runge-Kutta schemes for the large time step integration, and sound and gravity waves are treated implicitly in the vertical direction and explicitly in the horizontal directions for the small time step integration. Energy correction is introduced to ensure the conservation of total energy at every small time integration. The hydrological process including the warm rain cloud processes is also incorporated with the flux form equations. The scheme is tested with the non-hydrostatic core of the new global model.

Besides the conservations, two improvements have been devised in the model. First, we use more accurate formulas of the thermodynamics of the moist atmosphere by taking account of the effects of specific heats of water substances and the dependence of latent heat on temperature. These effects are generally neglected in most numerical models. Second, we introduced a conservative semi-Lagrangian method for the transportations due to rain.

We have performed experiments of a squall-line and direct calculation of a radiative-convective equilibrium to confirm the conservations of mass and energy. We found that, if the accurate moist thermodynamics are used, the total rain is reduced more than 10\% in the squall-line experiment in comparison to the case when the usual simplified thermodynamics are used. We also found that the change in energy due to transportation of rain are very large and cannot be negligible in the flux-form formulation, while that in momentum could be
negligible.

M. Shoucri, Institut de Recherche Hydro-Québec (IREQ),
A. Qaddouri, M.Tanguay and
J. Coté, Meteorological Service of Canada
A fractional steps method for the numerical solution of the shallow water equations

It has been pointed out some time ago by Yakimiw and Robert (Mon. Wea. Rev., 1986) that the method of fractional steps for the numerical solution of the shallow water equations has the advantage of reducing the multidimensional matrix inversion problem into an equivalent one-dimensional problem, so the technique becomes very simple and very attractive to apply provided it is accurate and stable enough.

We apply this fractional steps technique by splitting the shallow water equations, and successively integrating in every direction along the characteristics using the Riemann invariants (which are constant quantities along the characteristics), associated with cubic spline interpolation (Shoucri, J. Comp. Phys. ,1986). The method is tested on simple advection models as well as the full shallow-water equations, and shown to be accurate and stable. The linear analysis of the equations show the method is unconditionally stable, reproducing exactly the frequency of the slow mode, while the frequency of the fast modes is exact to second order.

Extension to the three-dimensional weather forecast equations will be discussed.

Reiji Suda, Nagoya University
Fast spherical harmonics transform of FLTSS and its evaluation

We are developing a library routine set called FLTSS (Fast Legendre Transform with Stable Sampling), with which the spherical harmonics transform can be done in time $O(T^2 \log T)$. This is the first report of developments and applications of the FLTSS.

Our algorithm is based on approximate fast polynomial interpolation using the FMM (Fast Multipole Method). FLTSS provides routines for fast evaluations and expansions of the associated Legendre functions and those derivatives. The evaluation points can be the Gauss points.

We will discuss two topics. The first topic is the performance. The speed-up rate in the floating operation counts reaches 9.1 for $T = 4095$ and $\epsilon = 10^{-6}$. However, our first implementation does not provide such high speed in CPU time, perhaps because of the fine granularity of computations. We are now developing several versions of the code with various schemes of tuning, and the results will be shown in the presentation.

The second topic is the error. Because of the use of the FMM, our algorithm has trade-off between the computational costs and the approximation error. We are evaluating the effects of the errors of the transform on the solutions using shallow water equations and non-divergent flow equations. Our code controls the error with weights about $n$ (i.e. requires lower precision for higher wave numbers), and the effects of the weights on the results will be also discussed.

Hirofumi Tomita, Masaki Satoh, and Koji Goto, Frontier Research System for Global Change
Development of global non-hydrostatic model using icosahedral grid

When the horizontal resolution becomes high in atmospheric general circulation model, the spectral method may have some difficulties
for high performance computing: the Legendre transformation and extensive data movement between computer nodes on a massively parallel computer. On the other hand, when the simple latitude-longitude grid is used in the grid method, the pole problem, occurs.
The grid spacing near the poles becomes very small as the resolution becomes high.

In order to radically overcome the pole problem, it is needed to use another type of grid which is distributed as homogeneously as possible on the sphere. The Next Generation Model Research Group in Frontier Research System for Global Change ( FRSGC ) has started to develop a new model based on the icosahedral geodesic grid, which is one of the quasi-homogeneous grids. Using the shallow water equations, we have overcome some difficulties in the icosahedral grid configuration[1].

When a simulation with high resolution in the horizontal direction is performed, the equation in the vertical dynamics should be also reconsidered from the usually used hydrostatic equation to the non-hydrostatic one. In almost of the existing non-hydrostatic models,
the conservations of several important quantities are not completely satisfied. Since our purpose is on the global climate simulation,
the conservations of mass and total energy are required. Our research group has been developing a new non-hydrostatic scheme, in which
such properties are completely satisfied[2]. We apply this non-hydrostatic scheme to the global model based on the icosahedral grid configuration.

In this Workshop, we introduce the formulation of used scheme in our non-hydrostatic global model and show the first result.
we also summarize the computational performance of our model on Earth Simulator.

[1] H.Tomita, M.Tsugawa, M.Satoh, and K.Goto (2001),
Shallow Water Model on a Modified Icosahedral Geodesic Grid
by Using Spring Dynamics,''
J. Compt. Phys., vol.174, p.579

[2] M.Satoh (2002),
Conservative scheme for the compressible non-hydrostatic models
with the horizontally explicit and vertically implicit time integration
scheme,''
Mon. Wea. Rev., in press.

Motohiko Tsugawa, Yukio Tanaka and Seong Young Yoon, Frontier Research System for Global Change
Development of a Global Ocean Model on Quasi-Homogeneous Cubic Grid

We are planning to apply eddy-resolving global ocean model to climate simulations on the Earth Simulator, a massively parallel vector computer. The model will be with about 0.1 deg resolution and longitude-latitude grid model will be adopted. However, in longitude-latitude models, meridians converge around the poles and it causes inefficient computation. The computational efficiency is critical in ocean climate simulations because it takes more than thousand years to spin up the global ocean. In order to construct numerically efficient model without deterioration of physical performance, we are exploring both numerical schemes and computational techniques.

We are now trying to introduce a quasi-uniform grid. Currently, we are using quasi-homogeneous cubic grids of Purser and Rancic(1998). This kind of cubic grid has non-orthogonal coordinate and has 8 singular points on the sphere. However, by using B-grid and by employing an adequate definition of metric tensors, we found that the dynamical core of the cubic grid is simple and treats those drawbacks appropriately. In the presentation we will talk about the development of global ocean model in FRSGC. We will discuss finite difference schemes on the cubic grid, treatment of the singular points and will show preliminary results of the cubic grid global model.

Agathe Untch and Mariano Hortal (ECMWF)
Tests with cubic spline interpolation in the vertical advection of the ECMWF semi-Lagrangian model in connection with a cubic finite-element discretization in the vertical

Recently ECMWF has implemented in its operational semi-Lagrangian model a finite element discretization for the vertical based on cubic B-splines. This provided the infrastructure also for the use of cubic splines for the vertical interpolations in the semi-Lagrangian advection scheme.

Semi-Lagrangian schemes using cubic splines are known to be less diffusive and less dispersive then schemes using cubic Lagarange interpolation but can suffer from noise due to large over- and undershooting near sharp gradients. We will present results from our investigations into the benefit or otherwise of using cubic spline interpolation for the vertical advection of the dynamical variables (winds, temperature and humidity) and for tracers in the ECMWF model and discuss ways of alleviating the noise problem associated with spline interpolation.

Nigel Wood, Met Office, J. Coté, Meteorological Service of Canada, Andrew Staniforth
A framework for the analysis of physics-dynamics coupling strategies

Co-authors Andrew Staniforth and Jean Cote.

NWP and climate models are viewed as consisting of two distinct
modules: the adiabatic, inviscid equation set (the dynamics); and the parametrised diabatic and viscous forcings (the physics). Each module is generally considered in isolation. The inevitable coupling between them is often done in a seemingly ad hoc manner - any justification usually being argued on physical grounds rather than founded in numerical analysis.

Recent developments mean that the dynamical components of the models are generally very accurate and numerically stable while considerable effort continues to be expended on improving the accuracy, and also the scope, of the various physical parametrisations. However, the recent work of Caya, Laprise and Zwack (1998) [CLZ98] well illustrates that the way the parametrisations are coupled to the dynamics can significantly limit the
accuracy and numerical stability of the model as a whole. Additionally they show that the choice of coupling can impact dramatically on the accuracy of the steady state that the model can achieve (which is generally the most understood and proven aspect of the physical parametrisations). It seems clear therefore that if the community is to reap the full benefits of its investment in improved dynamics and physics modules, more attention must be paid to analysing, understanding and improving the strategy for the physics-dynamics coupling.

The complexity of the physics schemes and the non-linear nature of the coupled system makes this a non-trivial task. However, CLZ98 have shown that useful progress can be made by judiciously reducing the problem to its essence. They presented a highly simplified model, or paradigm, of a physics-dynamics coupling and used it to diagnose the source of a problem in their global model. We have extended their paradigm to allow for the effects of uniform advection and any number of temporally and spatially varying forcings. These extensions permit the study of rather more complex physics-dynamics coupling strategies (as well as numerical issues such as spurious semi-Lagrangian orographic resonance) but in an analytically tractable framework. This extended framework will be presented as will results of its application to four coupling strategies in the context of a given spatio-temporal forcing coupled with a diffusive process.

K.S. Yeh, NASA Goddard Space Flight Center, M. Fox-Rabinovitz and S.J. Lin
A Variable-Resolution Finite-Volume General Circulation Model

The NASA/NCAR finite-volume General Circulation Model has been generalized to variable resolution with multiple areas of interest for both regional and global applications. As its original uniform-resolution version, the variable-resolution version also conserves the air mass, absolute vorticity and thermodynamic energy both locally and globally. One unique advantage of implementing variable resolution with finite-volume dynamics is the prevention of computational noise by the intrinsic diffusion associated withthe monotonicity constraints. While the numerical solution also benefits from the scale selectivity of the intrinsic diffusion in high-resolution areas, the diffusion effect in low-resolution areas is, however, overly excessive when the grid is aggressively stretched. A so-called "semi-stretched grid" is thus introduced to mitigate the deterioration of the numerical solution in low-resolution areas, and to reduce the amplitude of dispersion due to the irregularity of resolution. The semi-stretched grid is also designed with multiple areas of interest which can be used to focus the resolution on multiple sources of forcing. Various experiments have been conducted with the variable-resolution model, and preliminary results are quite encouraging.