Verification of a Modified Nonhydrostatic Global Spectral Dynamical Core Based on the Dry-Mass Vertical Coordinate: Three-Dimensional Idealized Test Cases

1.   Introduction

We have been developing a modified nonhydrostatic (NH) global spectral dynamical core using a dry-mass vertical coordinate (Peng et al., 2019). Importantly, it inherently solves the so-called non-conservation problem of the dry air mass at the level of the continuous set of equations, which exists in the original Aladin-NH like models [e.g., nonhydrostatic Integrated Forecast System (IFS) at ECMWF; Wedi et al., 2009] and other global models with vertical coordinates based on the moist hydrostatic pressure (e.g., Community Atmosphere Model (CAM) in Community Earth System Model version 1.0 (CESM1.0); Neale et al., 2012].

A notable feature of this modified dynamical core different from the original Aladin-NH like dynamical models (Bénard et al., 2010) is the combined use of a dry-mass vertical coordinate and a new temperature variable defined as ${T_{\rm m}} = {p \mathord{\left/ {\vphantom {p {({\rho _{\rm d}}{R_{\rm d}})}}} \right. } {({\rho _{\rm d}}{R_{\rm d}})}}$, where $p$ is the full pressure (vapor plus dry air), ${\rho _{\rm d}}$ is the density for dry air, and ${R_{\rm d}}$ is the gas constant for dry air. Physically, the variable ${T_{\rm m}}$ is the temperature that dry air would have if it has the same pressure as the moist air while its density remains unchanged, which is inspired by the modified potential temperature employed in the Model for Prediction Across Scales-Atmosphere (MPAS-A; Skamarock et al., 2012). This modified temperature can be further expressed as ${T_{\rm m}} = T\left[ {1 + ({R_{\rm v}}/{R_{\rm d}}){q_{\rm v}}} \right] \approx T(1 + 1.61{q_{\rm v}})$ with $T$ the temperature, ${R_{\rm v}}$ the gas constant for water vapor, and ${q_{\rm v}}$ the mass mixing ratio (mass per mass of dry air) of water vapor, and thus, it is different from the virtual temperature ${T_{\rm v}}$. As stated in Peng et al. (2019), taking ${T_{\rm m}}$ as state variable ensures that the mass continuity equation is naturally expressed in terms of dry-mass conservation.

A comparison between this modified dynamical core and the original Aladin-NH like dynamical kernel is briefly summarized in Table 1. The main motivation of such modifications on the original one is to conserve the mass of dry air rather than moist air, as it is just the fact in a precipitating atmosphere. Essentially, assuming that the total mass of moist air is conserved is equivalent to neglecting the mass sink effect associated with precipitation flux, which should be problematic for heavily precipitating systems such as tropical cyclone (Lackmann and Yablonsky, 2004; Li et al., 2023) and extreme rainfall event (Yang et al., 2023).

Table  1.  A comparison between the modified dynamical core and the original Aladin-NH like dynamical kernel. All symbols in the table are listed in the Appendix

Download as CSV

Model aspect	Modified dynamical core	Aladin-NH like one
Prognostic variable	${\pi _{\rm {ds}}}$, u, ${q_j}\left( {j = {\rm{v,c,r} }, \cdots } \right)$ ${T_{\rm{m}}} = T\left[ {1 + \left( { { { {R_{\rm{v}}} } /{ {R_{\rm{d}}} } } } \right){q_{\rm{v}}} } \right]$ $\hat q = \ln \left( { {p / { {\pi _{\rm{d}}} } } } \right)$ $\begin{aligned} & dl = d + \chi \\[-3pt] & d = - g\left( { {p / { {m_{\rm{d} } }{R_{\rm{d} } }{T_{\rm{m} } } } } } \right){ {\partial w}/ {\partial {\eta _{\rm{d} } } } }\\[-3pt] & \chi = \left( { {p/{ {m_{\rm{d} } }{R_{\rm{d} } }{T_{\rm{m} } } } } } \right)\nabla \phi \cdot { {\partial { {\boldsymbol u} } }/ {\partial {\eta_{\rm{d} } } } }\\[-3pt] & {D_3} = \nabla \cdot { {\boldsymbol u} } + d + \chi \end{aligned}$	${\pi _{\rm{s}}}$, ${\boldsymbol{u} }$, ${q_j}\left( {j = {\rm{v,c,r} }, \cdots } \right)$ ${T_{}}$ $\hat q = \ln \left( {{p/\pi }} \right)$ $\begin{aligned} & dl = d + \chi \\[-3pt]& d = - g\left( { {p /{m{R_{\rm{d} } }T} } } \right){ {\partial w} /{\partial \eta } }\\[-3pt]& \chi = \left( { {p/{mRT} } } \right)\nabla \phi \cdot { {\partial {\boldsymbol{u} } }/{\partial \eta } }\\[-3pt]& {D_3} = \nabla \cdot {\boldsymbol{u} } + \left( { { { {R_{\rm{d}}} } / R} } \right)d + \chi \end{aligned}$
Vertical coordinate	${\pi _{\rm{d}}} = A\left( { {\eta _{\rm{d}}} } \right) + B\left( { {\eta _{\rm{d}}} } \right){\pi _{\rm ds} }$ ${m_{\rm{d}}} = { {\partial {\pi _{\rm{d}}} } / {\partial {\eta _{\rm{d}}} } }, \;$${ {\partial {\pi _{\rm{d}}} }/{\partial z} } = - {\rho _{\rm{d}}}g$	$\pi = A\left( \eta \right) + B\left( \eta \right){\pi _{\rm{s}}}$ $m = {{\partial \pi } /{\partial \eta }} , \;$${{\partial \pi }/{\partial z}} = - \rho g$
Prognostic equation	$\dfrac{ { {\text{d} }{\boldsymbol{u} } } }{ { {\text{d} }t} } + \dfrac{ {\gamma {R_{\rm{d} } }{T_{\rm{m} } } } }{p}\nabla p + \dfrac{\gamma }{ { {m_{\rm{d} } } } }\dfrac{ {\partial p} }{ {\partial {\eta _{\rm{d} } } } }\nabla \phi = {F_{\boldsymbol{u} } }$ $\begin{gathered} \dfrac{ { {\text{d} }dl} }{ { {\text{d} }t} } + {g^2}\dfrac{p}{ { {m_{\rm{d} } }{R_{\rm{d} } }{T_{\rm{m} } } } }\dfrac{\partial }{ {\partial {\eta _{\rm{d}}} } }\left[ {\dfrac{\gamma }{ { {m_{\rm{d} } } } }\dfrac{ {\partial p} }{ {\partial {\eta _{\rm{d} } } } } - 1} \right] \\[-3pt] - g\dfrac{p}{ { {m_{\rm{d} } }{R_{\rm{d} } }{T_{\rm{m} } } } }\dfrac{ {\partial {\boldsymbol{u} } } }{ {\partial {\eta _{\rm{d} } } } } \cdot \nabla w - d\left( {\nabla \cdot {\boldsymbol{u} } - {D_3} } \right) \\[-3pt] = - g\dfrac{p}{ { {m_{\rm{d} } }{R_{\rm{d} } }{T_{\rm{m} } } } }\dfrac{ {\partial {F_{w} } } }{ {\partial {\eta _{\rm{d} } } } } + \dot \chi \\[-3pt] \end{gathered}$ $\dfrac{ { {\text{d} }{T_{\rm{m}}} } }{ { {\text{d} }t} } + \dfrac{ {R{T_{\rm{m}}} } }{ { {c_{{v}}} } }{D_3} = \dfrac{ { {c_{{p}}} } }{ { {c_{{v}}} } }{H_{\rm{m}}}$ $\dfrac{ { {\text{d} }\hat q} }{ { {\text{d} }t} } + \dfrac{ { {c_{ {p} } } } }{ { {c_{v} } } }{D_3} + \dfrac{ { { {\dot \pi }_{\rm{d} } } } }{ { {\pi _{\rm{d} } } } } = \dfrac{ { {c_{ {p} } } } }{ { {c_{ {v} } } } }\dfrac{ { {H_{\rm{m} } } } }{ { {T_{\rm{m} } } } }$ $\dfrac{ {\partial {\pi _{\rm ds} } } }{ {\partial t} } + \int_0^1 {\nabla \cdot \left( { {m_{\rm{d} } }{\boldsymbol{u} } } \right)} {\text{d} }{\eta _{\rm{d} } }{\text{ = 0} }$ $\dfrac{ { {\text{d} }{q_j} } }{ { {\text{d} }t} } = {F_{ {q_j} } };{q_j} = {q_{\rm{v}}},{q_{\rm{c}}},{q_{\rm{r}}}, \cdots$	$\dfrac{ { {\text{d} }{\boldsymbol{u} } } }{ { {\text{d} }t} } + \dfrac{ {RT} }{p}\nabla p + \dfrac{1}{m}\dfrac{ {\partial p} }{ {\partial \eta } }\nabla \phi = {F_{\boldsymbol{u} } }$ $\begin{gathered} \dfrac{ { {\text{d} }dl} }{ { {\text{d} }t} } + {g^2}\dfrac{p}{ {m{R_{\rm{d} } }T} }\dfrac{\partial }{ {\partial \eta } }\left[ {\dfrac{1}{m}\dfrac{ {\partial \left( {p - \pi } \right)} }{ {\partial \eta } } } \right] \\[-3pt] - g\dfrac{p}{ {m{R_{\rm{d} } }T} }\dfrac{ {\partial {\boldsymbol{u} } } }{ {\partial \eta } } \cdot \nabla w - d\left( {\nabla \cdot {\boldsymbol{u} } - {D_3} } \right) \\[-3pt] = - g\dfrac{p}{ {m{R_{\rm{d} } }T} }\dfrac{ {\partial {F_{{w} } } } }{ {\partial \eta } } + \dot \chi \\[-3pt] \end{gathered}$ $\dfrac{ { {\text{d} }T} }{ { {\text{d} }t} } + \dfrac{ {RT} }{ { {c_{{v}}} } }{D_3} = \dfrac{Q}{ { {c_{{v}}} } }$ $\dfrac{ { {\text{d} }\hat q} }{ { {\text{d} }t} } + \dfrac{ { {c_{{p}}} } }{ { {c_{{v}}} } }{D_3} + \dfrac{ {\dot \pi } }{\pi } = \dfrac{Q}{ { {c_{{v}}}T} }$ $\dfrac{ {\partial {\pi _{\rm{{\rm{s}}} } } } }{ {\partial t} } + \int_0^1 {\nabla \cdot \left( {m{\boldsymbol{u} } } \right)} {\text{d} }\eta {\text{ = 0} }$ $\dfrac{ { {\text{d} }{q_j} } }{ { {\text{d} }t} } = {F_{ {q_j} } };{q_j} = {q_{\rm{v}}},{q_{\rm{c}}},{q_{\rm{r}}}, \cdots$

 | Show Table

DownLoad: CSV

In this paper, we focus on the performance of this modified dynamical core in applications on the sphere by conducting three-dimensional (3D) benchmark cases as those defined in the context of the Dynamical Core Model Intercomparison Project (DCMIP; Ullrich et al., 2012, 2016, hereinafter DCMIP2012 and DCMIP2016), the reduced-radius sphere applications (Klemp et al., 2015), and the climate-focused evaluation methods (Held and Suarez, 1994; Thatcher and Jablonowski, 2016, hereinafter TJ2016). We firstly validate the NH global spectral dynamical core using some dry benchmark tests, including steady state, dry baroclinic waves, mountain waves in non-sheared and sheared background flows, and a dry Held–Suarez test. Without moisture, the results simulated by the modified and original Aladin-NH like ones should be identical. In order to confirm its capability on address the NH effects, a comparison with the results from the hydrostatic one is also made. Then, we compare the performance of the modified dynamical core and the original Aladin-NH one by using moist test cases with simple physics. As the splitting supercell test had been presented in Peng et al. (2019), only the moist baroclinic wave, idealized tropical cyclone, and moist Held–Suarez test cases are considered here.

The remainder of the paper is organized as follows. As a supplement to Section 2.3 of Peng et al. (2019), some more details on the temporal and spatial discretization are presented in Section 2. The simulation results for five dry test cases are presented in Section 3 and for three moist test cases in Section 4, respectively. Finally, a concluding summary is given in Section 5.

2.   On solving the govern equations

2.1   Spatial and temporal discretization

In our nonhydrostatic dynamical cores presented here, spatial discretization is implemented by a spectral method with reduced, linear Gaussian grid in the horizontal direction and a finite-difference scheme in the vertical direction, respectively. In order to largely alleviate its performance bottleneck, the spectral method employed here is a fast spherical harmonic transform algorithm developed by Yin et al. (2018, 2019), which is actually based on the pioneering research work of Tygert (2010). The basic principles for the vertical discretization follow Bubnová et al. (1995), which ensures that the vertical scheme is energy conserving and hydrostatic angular-momentum conserving (Bénard and Mašek, 2013). The vertical domain is divided into L layers labeled $(1, \cdots ,L)$ from top to bottom. Each layer is separated from its neighbors by interfaces labeled $(\tilde 0, \cdots ,\tilde L)$. The temporal discretization is implemented by the two-time-level semi-implicit semi-Lagrangian (2-TL SI SL) scheme described in Bénard et al. (2010), which allows a relatively large time step and therefore is more efficient.

2.2   The coupling to the physics

Some physical processes are necessary for the following test cases. The dynamical core is directly coupled with the simple-physics package from https://github.com/ClimateGlobalChange/DCMIP2016, including Kessler-type cloud microphysics, surface fluxes, and boundary-layer mixing. Moreover, the physics–dynamics coupling employs the SLAVEPP (Semi-Lagrangian Averaging of Physical Parameterizations) method described by Wedi (1999).

Additional spatial dissipations are also implemented in the dynamical core, including a combination of an implicit fourth-order horizontal diffusion (typically used in climate model configuration) in spectral space and the explicit horizontal second-order diffusion formulation of Smagorinsky (1963) in grid space, which is consistent with the standard MPAS dissipation configuration (Skamarock et al., 2012).

A detailed description about the implicit horizontal diffusion of order 2r in spectral space has been presented in Peng et al. (2019) [see also ECMWF (2016)], where 2r is the power of the horizontal diffusion scheme. For any model variable $\psi$, the constant diffusion coefficient is defined as:

with ${N_{\max }}$ the maximum truncation wavenumber and $a$ the earth radius. The operator order r is set to 2, the timescale is set as $\tau = 6\Delta t$, and the enhancing coefficient is set as ${h_\psi } = 1$ for relative vorticity, horizontal divergence, and vertical divergence, while ${h_\psi } = 5$ for pressure departure, temperature, and moist species.

The horizontal second-order diffusion formulation of Skamarock (1963) on model coordinate surfaces is expressed as

with an eddy viscosity K_h is determined from the horizontal deformation using a Smagorinsky first-order closure approach:

where the continuous deformation tensor components are defined as:

where $\phi = gz$ is the geopotential, and ${\phi _x} = {{\partial \phi } \mathord{\left/ {\vphantom {{\partial \phi } {\partial x}}} \right. } {\partial x}}$ and ${\phi _y} = {{\partial \phi } \mathord{\left/ {\vphantom {{\partial \phi } {\partial y}}} \right. } {\partial y}}$ are the zonal and meridian gradients of the geopotential, respectively. In Eq. (3), the length scale $l = {{2\pi a} \mathord{\left/ {\vphantom {{2\pi a} {\left( {2{N_{\max }} + 1} \right)}}} \right. } {\left( {2{N_{\max }} + 1} \right)}}$ for linear grids, ${C_{\rm{s}}}$ is a constant and in the following test cases is set to 0.125, which is consistent with the MPAS configuration for DCMIP2016 test (Skamarock et al., 2016). For scalar mixing, the eddy viscosity is multiplied by the inverse turbulent Prandtl number $P_r^{ - 1} = 3$.

3.   Dry test cases

In this section, we validate the nonhydrostatic global spectral dynamical core using some dry benchmark tests without any additional physics. Specifically, the steady state test is carried out to check the computational accuracy of the numerical approach employed by the dynamical core. The dry baroclinic wave test is used to check the performance of the dynamical core to a controlled, evolving instability. The mountain wave tests on a reduced-radius sphere are carried out to examine the nonhydrostatic responses of the dynamical core to flows over a mountain profile with/without linear vertical wind shear. The dry Held–Suarez test is conducted to evaluate the performance of the dynamical core in simulating statistical properties of the atmospheric general circulation.

3.1   Steady state

The underlying idea behind this test is to analyze how well the dynamical core maintains a steady, balanced initial state, which is an analytic solution to the hydrostatic primitive equation. The initial conditions for this test case are constructed as described in Jablonowski and Williamson (2006, hereinafter JW2006).

The horizontal resolution employed in the control experiment for the steady state test is T_L159, corresponding to approximate grid spacing of 125 km near the equator. There are 91 layers in the vertical with the last full model level at 0.01 hPa (~75 km), which are consistent with previously operational set of IFS (https://confluence.ecmwf.int/display/UDOC/Model+level+definitions, the same below). After the initialization, the simulations are run for 30 days, with fields output every day. To investigate the stability of the nonhydrostatic dynamical core and meanwhile explore the effects of time step and horizontal resolution, simulations with varying time steps from 600 to 3600 s at fixed resolution T_L159L91 and with varying horizontal resolutions from T_L39 to T_L639 at fixed 91 layers and time step of 1200 s are also conducted. Neither sponge layer nor any of additional spatial dissipations mentioned above is used.

The vertical cross-sections of temperature, potential temperature, zonal velocity, and relative vorticity on day 30 for the simulation at the resolution T_L159L91 with time step of 600 s are shown in Fig. 1. It is clear that there are nearly no visible differences between the simulated fields and the corresponding initial fields, that is to say, the present nonhydrostatic dynamical core can perfectly maintain the initial balanced state. We also analyze the other simulations with different time steps or with different horizontal resolutions (omitted). The corresponding results also show that each simulated state on day 30 is also nearly identical to the initial balanced state without any visible differences, even for the simulations with very large time step of 3600 s or with very coarse horizontal resolution T_L39.

Download

Larger image

PowerPoint slide

Fig  1.  Vertical cross-sections of (a) temperature, (b) potential temperature, (c) zonal velocity, and (d) relative vorticity on day 30 for the simulation at the resolution T_L159L91 with time step of 600 s in the steady state test.

To further quantify the extent to which the simulated results deviate from the analytic solution, we also calculate the ${l_{\text{2}}}$ (root-mean-square) error norm of the zonal-wind field $u$. The expression of the corresponding ${l_{\text{2}}}$ error norm is given by

where the overbar $\overline {\left( {} \right)}$ represents the zonal average, ${w_j}$ is the Gaussian weight at the latitude points ${\varphi _j}$, and $\Delta {\eta _k} = {\eta _{k + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} - {\eta _{k - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}$ denotes the thickness of a model layer. The results calculated by using Eq. (5) for the different simulations with varying time steps and with varying horizontal resolutions are presented in Figs. 2a, b. For all the simulations, the ${l_2}$ error norms grow slightly after the start of simulation, which is mainly due to the initial adjustment of the zonal mean by the nonhydrostatic dynamical core. Thereafter, the ${l_2}$ errors exhibit nearly flat profiles with almost no increase over time. The maximum of the ${l_2}$ errors is almost no more than 0.005 m s⁻¹, which is only a quarter of that simulated with the Eulerian (EUL) dynamical core by JW2006 (see their Fig. 4a). That is to say, the convergence of ${l_2}$ errors as horizontal resolutions increase (at fixed time step) and as time steps decrease (at fixed horizontal resolution) is very strong. It is concluded that the nonhydrostatic dynamical core can maintain the balanced steady state very well. These results also validate the computational accuracy of the numerical approach employed by the dynamical core.

Download

Larger image

PowerPoint slide

Fig  2.  Time series of root-mean-square ${l_2}$ norm of $\bar u - {\bar u_{t = 0}}$ for the different simulations with (a) varying time steps at the resolution T_L159L91 and (b) varying horizontal resolutions and 91 levels at time step of 1200 s in the steady state test.

3.2   Dry baroclinic waves

Following previous studies (JW2006; Ullrich et al., 2012), the evolution of baroclinic waves is triggered by overlaying the above steady-state initial conditions with the Gaussian-type zonal wind perturbation in the Northern Hemisphere, expressed as

where ${R_p} = {a \mathord{\left/ {\vphantom {a {10}}} \right. } {10}}$ is the radius of perturbation, $a$ is the earth radius, ${u_p} = 1{\text{ m }}{{\text{s}}^{ - 1}}$ is the maximum amplitude, and $r$ is the great circle distance relative to the perturbation center $\left( {{\lambda _{\rm{c}},{\varphi _{\rm{c}}}}} \right) = \left( {{{\text{π }} \mathord{\left/ {\vphantom {{\text{π }} 9}} \right. } 9},{{2{\text{π }}} \mathord{\left/ {\vphantom {{2{\text{π }}} 9}} \right. } 9}} \right)$.

After the initialization, the simulation of baroclinic waves is run at the resolution T_L319L91 for 10 days, with output every 3 hours. The time step employed here is $\Delta t = {\text{1200 s}}$, and the other configurations are identical to the steady state test case.

Figure 3 shows the evolutions of 850-hPa temperature and surface pressure in the simulated baroclinic wave life cycle. During the first 4 days, the wave grows very slowly. On day 6, two small-amplitude but visible waves appear in the temperature field, and the surface pressure shows two weak low- and high-pressure systems at almost the same places, propagating eastward. On day 8, the two high- and low-pressure systems have significantly deepened (the scales of systems also have become larger) and the waves in the temperature field have almost peaked, with the cold-air mass on the western side beginning to invade into the warm wave ridges. In addition, the third upstream wave has become visible. The wave breaking onsets on day 10, and three closed cells with relatively higher temperature form.

Download

Larger image

PowerPoint slide

Fig  3.  Evolution of the baroclinic waves from days 4 to 10 for (a, c, e, g) temperature (K) at 850 hPa and (b, d, f, h) surface pressure (hPa) simulated with the modified nonhydrostatic dynamical core at the resolution T_L319L91 and the time step ∆t = 1200 s.

Download

Larger image

PowerPoint slide

Fig  4.  As in Fig. 3, but for the 850-hPa relative vorticity (10⁻⁵ s⁻¹) on (a) days 7 and (b) 9.

The evolution of baroclinic wave life cycle presented here is almost identical to that reported in JW2006 (see their Fig. 5). More specifically, the low surface pressure centers simulated here on days 8 and 10 are as strong as those by the finite-volume (FV) dynamical core at the horizontal resolution 0.5° × 0.625° in JW2006 (see their Figs. 5c, d) and the fronts in the temperature field simulated here are as sharp as those in the latter correspondingly (see their Figs. 5g, h). Such similarities are consistent with the comparable horizontal resolution between these two. In addition, also analyzed are the simulated results in the Southern Hemisphere (figure omitted), indicating that there are no obvious grid-imprinting errors.

Download

Larger image

PowerPoint slide

Fig  5.  Convergences of the modified nonhydrostatic dynamical core with increasing horizontal resolution and 91 levels at fixed time step ∆t = 1200 s for (a, c, e, g) temperature (K) at 850 hPa and (b, d, f, h) surface pressure (hPa) on day 9.

Figure 4 displays the simulated 850-hPa relative-vorticity fields on days 7 and 9 by the modified nonhydrostatic dynamical core. Once again, the overall vorticity patterns simulated here agree very well with those reported for the FV dynamical core in JW2006 (see their Fig. 8c), even in the vicinity of the rolled up tongues of the vortices on day 9, although the former looks slightly stronger at some vorticity centers.

To further evaluate numerical convergence of the simulated baroclinic wave by the present dynamical core, a series of sensitivity experiments at varying horizontal resolutions are conducted. Figure 5 compares the baroclinic solutions on day 9 for different horizontal resolutions. The T_L159 run does not capture the sharpness of the fronts in the temperature field and the strength of the closed cells in the surface pressure as in the higher-resolution simulations. The T_L319 run is visually almost identical to the T_L639 simulation, although differences at the smallest scales can be seen as expected.

As an extremely important issue, the parallel performance of the modified nonhydrostatic dynamical core is briefly evaluated with dry baroclinic wave simulations. The run times at the second-highest and highest horizontal resolutions (T_L319 and T_L639, respectively) and 91 levels are listed in Table 2. The run time data are measured by the wallclock time needed to complete one model day on a 64-processor node of an Intel Xeon architecture. Note that only a pure MPI (Message Passing Interface) parallelization approach is considered here. Compared to the T_L319L91 simulation with time step 1200 s, the T_L639L91 simulation with time step 600 s corresponds to an increase in the workload by a factor of 8. It is shown that the increase in the wallclock time for the T_L639L91 simulation lies around a factor of 6.35, which implies that the modified nonhydrostatic core shows a superlinear parallel speedup as the horizontal resolution is doubled and the time step is halved.

Table  2.  Wallclock time (s) for one model day measured on a 64-processor node of an Intel Xeon architecture

Download as CSV

Horizontal resolution	Time step	Wallclock time for one day
TL319	1200 s	200 s
TL639	600 s	1269 s
Increase in run time by a factor of		6.35

 | Show Table

DownLoad: CSV

3.3   Mountain waves in non-sheared background flow

The mountain wave test cases proposed by Klemp et al. (2015) are conducted to evaluate the performance of the nonhydrostatic dynamical core in handling complex topography. These test cases, referred to as variations from the original cases in DCMIP2012, simulate the flow of an atmosphere with constant wind and stability over a specified mountain profile on a non-rotating, reduced-radius sphere. We also consider two different mountain profiles, a quasi-2D ridge, and a circular mountain.

Following the DCMIP2012 test case 2.1 [Eq. (7)], the formulation of the circular mountain profile is expressed as:

where $r$ is also the great circle distance but relative to the mountain center $\left( {{\lambda _{\rm c}},{\varphi _{\rm c}}} \right) = \left( {{{\text{π }} \mathord{\left/ {\vphantom {{\text{π }} 4}} \right. } 4},0} \right)$, defined as

Here, $X$ is the reduced-radius earth scaling factor, ${h_0} = 250{\text{ m}}$, ${d_0} = 5000{\text{ m}}$, and ${\xi _0} = 4000{\text{ m}}$.

In order to make the 3D test case similar to the 2D test case presented by Schär et al. (2002), Klemp et al. (2015) modified the circular mountain profile to a quasi-2D ridge profile, expressed as

with ${r_0}\left( \lambda \right) = \left( {\lambda - {\lambda _{\rm c}}} \right){a \mathord{\left/ {\vphantom {a X}} \right. } X}$. With this ridge-like terrain, the mountain-wave structure along its centerline (equator) should be almost identical to the 2D solutions presented in Schär et al. (2002).

The initial state of the atmosphere has a zonal wind in solid body rotation with no vertical wind shear, given by

with ${u_{{\rm {eq}}}} = 20{\text{ m }}{{\text{s}}^{{- 1}}}$. As pointed out by Klemp et al. (2015), in absence of vertical wind shear, the initial state of the atmosphere should be isothermal, i.e., $T\left( {\lambda ,\varphi ,z} \right) =$${T_{\rm {eq}}}$. In present test cases, ${T_{{\rm {eq}}}}$ is set to 300 K.

Using the hydrostatic and gradient wind balance relationships, we can further obtain initial pressure field, given by

The desired surface pressure ${p_{\rm s}}\left( {\lambda ,\varphi } \right)$ is obtained by setting $z = {z_{\rm s}}$ in the above equation.

After the initialization, the simulations with different mountain profiles are configured accordingly and integrated for 2 h. The reduced-radius scaling factor is set to $X = 166.7$, which minimizes the curvature effects and prevents mountain waves from circling the globe during a 2-h integration. There are also 91 layers in the vertical with the highest full model level at 0.01 hPa (~100 km), which is much higher than the model top used in previous studies (e.g., Schär et al., 2002; Klemp et al., 2015). Thus, no explicit absorbing layer is employed here. In order to suppress subgrid-scale noises, we configure present simulations with only the implicit fourth-order horizontal diffusion in spectral space, which is consistent with the fourth-order hyperdiffusion associated with the third-order transport scheme used in Klemp et al. (2015). Specially, in order to highlight the nonhydrostatic effects represented by the dynamical core, the simulations with nonhydrostatic and hydrostatic configurations for each mountain profile are conducted.

For the quasi-2D mountain ridge case, we first make a comparison between the nonhydrostatic and hydrostatic simulations at the resolution T_L159L91 and time step ∆t = 12 s. This choice leads to a grid spacing of approximate 1.1° (~720 m along the equator), which is consistent with Klemp et al. (2015) for the MPAS simulation. The vertical cross-sections of vertical velocity along the ridge centerline from these T_L159L91 simulations for the ridge-like terrain with h₀ = 250 m at t = 2 h are shown in Figs. 6a, b, respectively. For reference, the corresponding 2D linear analytic solution calculated as in the Appendix A of Klemp et al. (2015) is also shown in Fig. 6d. One can see that the nonhydrostatic core does simulate weak-amplitude mountain waves over the leading terrain profile (Fig. 6a), which is just the case represented by the 2D linear analytic solution (Fig. 6d); while the hydrostatic core simulates much stronger disturbances with very different structures (Fig. 6b). However, at the resolution T_L159L91, there is still a significant distortion of the steady mountain wave solved by the nonhydrostatic core. As discussed by Klemp et al. (2003), this distortion should be mainly caused by the numerical inconsistency in treating the metric terms associated with the terrain-following coordinate transformation. In present nonhydrostatic core, the advection operator is a semi-Largrangian one, with a third-order interpolation of variables to the departure points; while the metric terms associated with the vertical coordinate transformation are treated through spectral discretization, with an obscure order of accuracy depending on the maximum truncated wavenumber. Further comparison with Fig. 3 of Klemp at al. (2015) shows that the distortion evident in Fig. 6a here is very close to that in their Fig. 3c using fourth-order advection together with second-order numerics for the metric term, implying that this error may be mainly attributed to insufficient order of accuracy of the metric terms. In view of this, we then reperformed the simulation at higher spatial resolution of T_L319L91. As expected, the corresponding solution depicted in Fig. 6c, with higher resolution, largely reduces the distortion and thus better reproduces the analytic solution shown in Fig. 6d. Therefore, for the ridge terrain profile, it is spectral resolution T_L319 (~360 m) rather than T_L159 (~720 m) that is sufficient for the present nonhydrostatic core to accurately capture the wave development of interest above the mountain. In addition, it should also be noted that there exist some trailing perturbations between the nonhydrostatic results and the linear analytic solution at higher levels, which is due to spherical effects and is also observed from the MPAS simulation in Fig. 3 of Klemp at al. (2015).

Download

Larger image

PowerPoint slide

Fig  6.  Vertical cross-sections of vertical velocity (m s⁻¹) along the equator (centerline) from the (a) nonhydrostatic (nh) and (b) hydrostatic (hy) simulations at the resolution T_L159L91 and time step ∆t = 12 s for the ridge-like terrain with a maximum height h₀ = 250 m at t = 2 h; (c) as in (a), but for the nonhydrostatic simulation at the higher resolution T_L319L91; and (d) vertical cross-section of vertical velocity (m s⁻¹) from the 2D linear analytic solution for h₀ = 250 m [cf. Appendix A of Klemp et al. (2015)]. To facilitate comparison, the units of the x-coordinate have been converted from meter to degree in the linear analytic solution, neglecting the effects of spherical curvature.

For the circular mountain case, there are more significant smaller-scale structure in the nonhydrostatic wave train, as documented by Klemp et al. (2015). Therefore, the relatively higher horizontal resolution and smaller time step are suggested for this test case. Maintaining similarity with Klemp et al. (2015), we specify a higher spatial resolution T_L319L91 and a smaller time step ∆t = 6 s. The corresponding simulated results for the circular-mountain terrain are shown in Fig. 7. For comparison, we directly refer to the 3D linear analytic solution for the circular-mountain terrain published in Klemp et al. (2015, their Fig. 5f). It is also clearly found that only the nonhydrostatic core can capture the mountain-wave structure fairly well, and the results resolved by the present nonhydrostatic core are quite similar with the numerical results simulated by MPAS and even closer to the corresponding 3D linear analytic solution in Klemp et al. (2015).

Download

Larger image

PowerPoint slide

Fig  7.  As in Figs. 6a, b, but for the circular-mountain terrain at the resolution T_L319L91 and time step ∆t = 6 s.

3.4   Mountain waves in sheared background flow

Simulating mountain waves in sheared background flow originally introduced by Wurtele et al. (1987) serves as another particularly discriminating tool to evaluate the performance of the dynamical core (Keller, 1994; Wedi and Smolarkiewicz, 2009). In contrast to the tilted, vertically propagating mountain gravity wave above, the correct solution for this case is that of trapped, horizontally propagating gravity waves. We initialize the present case as in Wedi and Smolarkiewicz (2009) except that the location of the mountain center is shifted by 180 degrees in longitude. The mountain profile is of 3D elliptic shape, defined by

with ${l_\lambda } = \dfrac{a}{X}\arccos \left[ {{{\sin }^2}{\varphi _{\rm c}} + {{\cos }^2}{\varphi _{\rm c}}\cos \left( {\lambda - {\lambda _{\rm c}}} \right)} \right]$ and ${l_\varphi } =$$\dfrac{a}{X}\arccos \left[ {\sin {\varphi _{\rm c}}\sin \varphi + \cos {\varphi _{\rm c}}\cos \varphi } \right]$. Here, ${L_\lambda } = 2.{\text{5 km}}$ is the mountain half-width, ${L_\varphi } = {\left| {L_\lambda ^2 - L_f^2} \right|^{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}$denotes the meridional extent of the elliptic, and ${L_f}$ is the focus point distance with the mountain center $\left( {{\lambda _{\rm c}},{\varphi _{\rm c}}} \right) = \left( {{\text{π} \mathord{\left/ {\vphantom {\pi 2}} \right. } 2},0} \right)$ and the focus point $\left( {{\lambda _{\rm d}},{\varphi _{\rm d}}} \right) = \left( {{\text{π} \mathord{\left/ {\vphantom {\text{π} 2}} \right. } 2},{\text{π} \mathord{\left/ {\vphantom {\text{π} 3}} \right. } 3}} \right)$.

The background state of atmosphere is characterized by constant stability with Brunt–Väisälä frequency $N = 0.01{\text{ }}{{\text{s}}^{{ - 1}}}$, and zonal wind with a linearly vertical shear profile $u\left( {\varphi ,z} \right) = {U_0}\left( {1 + cz} \right)\cos \varphi$ below the tropopause (about 10.5 km) and constant in the stratosphere, where ${U_0}$ and $c$ are set to $10{\text{ m }}{{\text{s}}^{{- 1}}}$ and $2.5 \times {10^{ - 4}}{{\text{ m}}^{- 1}}$, respectively. For simplifying the specification of a constant stability, the initial atmosphere is assumed to be isothermal with ${T_0} = {{{g^2}} \mathord{\left/ {\vphantom {{{g^2}} {\left( {{c_p}{N^2}} \right)}}} \right. } {\left( {{c_p}{N^2}} \right)}}$ (Wedi and Smolarkiewicz, 2009). It should be pointed out that although the isothermal initial condition is inconsistent with the vertically sheared zonal flow due to the gradient wind relationship, it is still okay to use the unbalance initial state for the present case.

The numerical setup is based on Wedi and Smolarkiewicz (2009). The simulations are configured with a T_L159 horizontal resolution and 91 vertical levels. The reduced-radius scaling factor is set to $X = 321$, corresponding to a reduced radius ~20 km. In order to ensure approximately uniform vertical spacing in the troposphere, we employ the exponential stretch vertical coordinate, which is equivalent to a hybrid terrain-following coordinate ${\eta _{\rm d}}$ with $A\left( {{\eta _{\rm d}}} \right) = 0$ and $B\left( {{\eta _{\rm d}}} \right)$ being expressed as

The control parameters in Eq. (13) are set as ${N_{\rm s}} = 290$ and ${a_{\rm s}} = 0.5$, resulting in a nearly constant vertical spacing of $\Delta z \approx 210{\text{ m}}$ below the height of 12 km, which are comparable to that used in Wedi and Smolarkiewicz (2009). In order to avoid gravity wave reflecting off the model top, Rayleigh damping is applied to the vertical velocity above the height of 20 km (Klemp et al., 2008). The temperature of the reference state for the semi-implicit scheme is set to 1000 K. After the initialization, the simulations with nonhydrostatic and hydrostatic configurations respectively are run for 2 h with time step $\Delta t = 5{\text{ s}}$. The corresponding results are shown in Fig. 8.

Download

Larger image

PowerPoint slide

Fig  8.  Vertical cross-sections of vertical velocity (m s⁻¹) along the equator from the (a) nonhydrostatic (nh) and (b) hydrostatic (hy) simulations at t = 2 h for linearly sheared flow past a 3D elliptic mountain on the reduced-radius with X = 321.

As shown in Fig. 8, the nonhydrostatic dynamical core can resolve the trapped, horizontally propagating gravity waves well. Specifically, there is a trail of closed cells behind the mountain with a roughly horizontal wave-length of $\dfrac{{60}}{{180}}{\text{π }} \dfrac{a}{X} \; \times \; \dfrac{ 2}{ 3} \approx 14 \; {\rm{km}}$, which exhibits good quantitative agreement with the numerical solution in their Fig. 1 of Wedi and Smolarkiewicz (2009) and with the linear analytic solution of a similar case in Wurtele et al. (1987). In addition, the amplitudes of the vertical velocity simulated by nonhydrostatic core are also comparable to the numerical results in previous studies (e.g., Wedi and Smolarkiewicz, 2009).

3.5   Dry Held–Suarez test

As a climate-focused evaluation method, the dry Held–Suarez test with highly idealized forcing and dissipation, namely, Newtonian relaxation of the temperature field and Rayleigh damping of the low-level winds (Held and Suarez, 1994), is conducted to evaluate the performance of the nonhydrostatic dynamical core in simulating statistical properties of the atmospheric general circulation.

The Held–Suarez simulation is started from an isothermal state at rest, with some small perturbations added into the vorticity field to break the symmetry. The model is then integrated for 1200 days at the resolution T_L159L31 and time step $\Delta t = {\text{1800 s}}$. The model fields are output every 6 hours, and statistics are calculated by averaging the last 1000 days of integration. The implicit fourth-order horizontal diffusion in spectral space is activated to suppress subgrid-scale noises. No explicit sponge layer is employed here. As there is no analytic solution to the Held–Suarez test, the simulated results are directly compared to the CAM-EUL results from the official website at https://www2.cesm.ucar.edu/models/simpler-models/held-suarez.html.

Figure 9 presents the 1000-day mean climatology for various fields from the dry Held–Suarez simulation with the modified nonhydrostatic dynamical core. One can see that the modified dynamical core almost perfectly reproduces the CAM-EUL results on a quadratic Gaussian grid at T85L30 resolution, including the position, structure, and strength of both the westerlies in the midlatitudes and the easterlies in the tropics and polar regions, the eddy temperature variance as well as the meridional eddy momentum and heat fluxes. In addition, there are some visible differences in the lower stratosphere above 100 hPa, which may be related to the absence of the absorption layer near the top of the model.

Download

Larger image

PowerPoint slide

Fig  9.  Zonal-mean outputs from the dry Held–Suarez test with the modified dynamical core at the resolution T_L159L31 and time step $\Delta t = {\text{1800 s}}$: (a) zonal wind, (b) eddy temperature variance, (c) northward eddy momentum flux, and (d) northward eddy heat flux. All are 1000-day mean from days 200 to 1200 of integration.

Apart from the above statistics, we further evaluate the performance of the modified dynamical core in mass conservation by calculating the vertically integrated quantities of full mass $m$, total column water mass ${m_{\rm q}}$, and dry mass ${m_{\rm d}}$ as follows (Berrisford et al., 2011; Peng et al., 2019):

where ${q_{\rm t}} = {q_{\rm v}} + {q_{\rm c}} + {q_{\rm r}} + \cdots$ is the total mass mixing ratio of moist species. Note that, for the dry Held–Suarez test, there is no moisture and thus only the dry mass actually needs to be analyzed. Figure 10 presents the time change of global mean dry mass from day 0 during the 1200-day integration. It is shown that when no mass fixer is applied, the dry mass increases linearly with time, which is consistent with the performance of the IFS simulation [Fig. 4 in Malardel et al. (2019)]. Such behaviour should be mainly attributed to the non-conservative numerical treatment (i.e., the semi-Lagrangian advection scheme; Wong et al., 2013), and needs further consideration in future work.

Download

Larger image

PowerPoint slide

Fig  10.  Time change of global mean dry-air mass from day 0 during the 1200-day integration simulated with the modified dynamical core. The mass is defined as vertical integration through a whole atmospheric column and the global mean is computed with the Gaussian weights.

4.   Moist test cases

In this section, three additional moist benchmark tests are conducted to further evaluate moist aspects of the nonhydrostatic global spectral dynamical core. Both the results simulated from the modified and original cores are provided and compared in order to explore the significant effects of our modifications of the dynamical core on the development of precipitating systems.

4.1   Moist baroclinic waves

Moist baroclinic wave test is proposed to study the moisture feedbacks on the dynamical core. The present case is initialized as described in DCMIP2016, where the initial reference state is an analytical solution of the motion equations in height coordinate with a flat topography and zero velocity at the model surface, and the zonal wind perturbation takes the form of a simple exponential bell with a vertical taper. After the initialization, the simulation of moist baroclinic waves is run for 20 days at the resolution T_L159L91, comparable in horizontal resolution to the grid used by the reference results obtained by the 1° CAM-FV simulation from the website at https://www2.cesm.ucar.edu/models/simpler-models/fkessler/index.html. Only the Kessler-type cloud microphysics from the simple-physics packages is taken into account in present simulations. In order to suppress possible strong convection overturning, a relatively smaller time step of $\Delta t = {\text{300 s}}$ is employed in this moist case. Other configurations are identical to those described in the above test of dry baroclinic waves. We first make a comparison between the surface pressure obtained by present simulations with the available reference using CAM-FV on day 10 (figure omitted), highlighting the close resemblance between them, and then focus on the results of a longer integration time.

Figure 11 presents the results for the simulated moist baroclinic waves on day 15. It is shown that the solutions obtained from the modified dynamical core and the original one show close agreement with slight differences, even in this much later stage of the nonlinear baroclinic instability evolution. The analysis of the simulations is supplemented by comparing the temporal evolution of minimum surface pressure. In Fig. 12, one can see that the minimum surface pressure simulated by these two dynamical cores is nearly identical during the first 17 days, and subsequently, the differences between them become slightly obvious (left column).

Download

Larger image

PowerPoint slide

Fig  11.  Temperature at 850 hPa (left column; K) and surface pressure (right column; hPa) in the simulated moist baroclinic waves on day 15: the first two rows show the results for the modified dynamical core and the original one, respectively; while the last row shows the corresponding difference between the solutions.

Download

Larger image

PowerPoint slide

Fig  12.  Moist-precipitating baroclinic instability: (a) time series of minimum surface pressure (hPa) simulated by the modified dynamical core (red line) and the original one (blue line) and (b) accumulated precipitation (mm) through day 15 for the modified dynamical core.

In conclusion, although including moisture, the baroclinic instability solutions obtained with these two dynamical cores show no significant differences, which is consistent with the fact that the particular test case focuses on the large-scale hydrostatic regime and the water mass loss due to the large-scale precipitation here is still relatively small (right column).

4.2   Idealized tropical cyclone

The idealized tropical cyclone (TC) test case is recommended as a moist, deterministic test case of intermediate complexity, and its purpose is to simulate a rapid intensification of TCs by initializing an analytic vortex in a background environment (Reed and Jablonowski, 2011, 2012). The present case is initialized as described in the DCMIP2016 TC test case 162-bryanpbl, where the modified Bryan PBL scheme was used. After the initialization, the simulation of idealized TC with the full simple-physics package and additional spatial dissipations is performed at T_L639 resolution (i.e., approximate grid spacing of 31 km near the equator) with 31 vertical levels for 10 days. The PBL scheme used here is also the modified Bryan one. Both the simulations by the modified and original dynamical cores are conducted. As mentioned in the introduction, the precipitation mass sink effect is naturally included in the modified dynamical core, thus it is expected to obtain simulated results with lower pressure, stronger wind speeds, and heavier precipitation (Lackmann and Yablonsky, 2004). Moreover, it is also expected to get results closer to those from MPAS model presented in DCMIP2016 due to their similarities in dynamical core [see Peng et al. (2019) for a detailed explanation] and the nearly identical configurations for TC simulation.

Figures 13 and 14 display the overall characteristics of the simulated TC on day 10 by these two dynamical cores in a manner similar to previous results for other models (DCMIP2016), respectively. It is shown that both the modified and original dynamical cores can simulate a storm resembles a TC with very intense wind speeds, a relatively calm eye, a warm-core, and extreme precipitation, consistent with previous results for other models in DCMIP2016. However, there is also obvious difference between the simulated TCs in cyclone intensity and location. For the modified dynamical core, the maximum wind speed at 1500 m is 72.55 m s⁻¹ and the maximum temperature anomaly at 5000 m is 10.45 K; while for the original one, they are 61.17 m s⁻¹ and 9.00 K, respectively. That is to say, the idealized TC on day 10 simulated by the modified dynamical core is much stronger than that by the original one indeed. Further comparing with the results from MPAS model presented in DCMIP2016, one can found that the modified dynamical core can simulate an idealized TC with comparable intensity exceeding 70 m s⁻¹ at a height of 1500 m; while the original one produces an obviously weak one.

Download

Larger image

PowerPoint slide

Fig  13.  Snapshots of the simulated idealized tropical cyclone on day 10 by the modified dynamical core: (a) horizontal wind speed at a height of 1500 m, (b) zonal velocity at a height of 100 m, (c) temperature anomaly at a height of 5000 m, and (d) accumulated precipitation through day 10.

Download

Larger image

PowerPoint slide

Fig  14.  As in Fig. 13, but for the original dynamical core.

The time evolution of the simulated TCs represented by the minimum surface pressure is compared in Fig. 15. It is shown that both the two dynamical cores can capture the rapid intensification of the vortex during the first 6 days (Fig. 15a), and the apparent north-northwestward movement of the TC (Fig. 15b), which is consistent with previous results for other nonhydrostatic global models. However, there is also obvious difference between the dynamical cores in both the intensity and location of cyclone. Specifically, at test termination of 10 days, the modified dynamical core produces a stronger TC with minimum surface pressure around 937 hPa, which is 20 hPa lower relative to the original core simulation and thus much closer to the MPAS simulations (i.e., 930 hPa for 60-km resolution, 920 hPa for 30-km resolution, and 930 hPa for 15-km resolution); and the storm simulated by the former moves northward faster than that by the latter, resulting in an offset of approximately 4 degrees by 10 days.

Download

Larger image

PowerPoint slide

Fig  15.  Time evolution of the (a) minimum surface pressure and (b) location of the tropical cyclone simulated by the modified dynamical core (red) and the original one (blue). The location of the tropical cyclone is tracked by the minimum surface pressure.

4.3   Moist Held–Suarez test

TJ2016 developed a moist Aqua-planet variant of the Held–Suarez test, which is proposed to shed light on the nonlinear dynamics–physics moisture feedbacks. The simplified moist physics parameterizations for this moist case follow those from the above short-term TC test case with several modifications. These modifications make the simplified moist physics package appropriate for long-term climate studies; details refer to TJ2016. The initialization of the moist Held–Suarez simulation is based on the shallow atmosphere version of the baroclinic wave test for dynamical cores developed by Ullrich et al. (2014), which provides a steady-state moist atmosphere without topography. The moist surface pressure is set to ${p_{\rm s}} = {p_0} = 1000{\text{ hPa}}$ everywhere, and the initial specific humidity profile is specified with the maximum specific humidity ${q_0} = 18{\text{ g k}}{{\text{g}}^{- 1}}$, the horizontal half-width with latitude ${\phi _{{\rm hw}}} = {{2{\text{π }}} \mathord{\left/ {\vphantom {{2{\text{π}}} 9}} \right. } 9}$, and the vertical half-width with pressure ${p_{{\rm hw}}} = 340{\text{ hPa}}$ (details refer to Appendix A of TJ2016). Other configurations for the moist Held–Suarez test are identical to the dry Held–Suarez test in Section 3. The zonal mean climatology is defined as the average of days 210–1200 of the simulations using the 6-h interval history fields.

Figure 16 compares the time-mean zonal-mean zonal wind, temperature, and precipitation rate simulated by the modified and original dynamical cores, respectively. Overall, the zonal-mean climatologies simulated by these two cores are quite similar and agree quite well with the climatology obtained by the CAM-FV simulation at approximately 2° resolution from the official website https://www2.cesm.ucar.edu/models/simpler-models/moist_hs/index.html. There are nearly no visible differences in the temperature fields simulated by these two cores (Figs. 16b, e). An interesting zonal-wind difference between these two simulations is the horizontal symmetry of the easterlies (see the −4-m s⁻¹ contour) in the polar regions (Figs. 16a, d), which maintains perfect in the simulation with the modified dynamical core. The most significant difference lies in the zonal mean precipitation rate in the tropics (Figs. 16c, f), which is obviously enhanced by the modified dynamical core. Specifically, the equatorial peak precipitation rate exceeds 25 mm day⁻¹ in the simulation with the modified dynamical core; while it does not exceed 21 mm day⁻¹ in the simulation with the original one.

Download

Larger image

PowerPoint slide

Fig  16.  Zonal-mean climatologies from the moist Held–Suarez simulations at the resolution T_L159L31 and time step $\Delta t = {\text{1800 s}}$ with (upper panels) the modified dynamical core and (lower panels) the original one: (a, d) zonal wind, (b, e) temperature, and (c, f) precipitation rate. All are 990-day means from days 210 to 1200 of integration.

The time-mean zonal-mean vertical velocities in the tropics are shown for the modified and original cores in Figs. 17a, c, respectively. Both simulations show a narrow updraft area close to the equator, divided into two branches above 200 hPa. This corresponds to the upward branch of the Hadley circulation. The original dynamical core anchors the peaks in the lower atmosphere near 800 hPa, which is consistent with the moist idealized simulations without the deep convection parameterization in TJ2016 (see their Figs. 2a, b). However, the modified dynamical core significantly enhances the updraft speeds, and extends the peaks further up. As a result, the vertical velocity shows two equatorial updraft peaks in the simulation by the modified one, the first one in the lower atmosphere near 800 hPa and the second one in the upper atmosphere near 400 hPa, which is largely consistent with the Aqua-planet simulation with the deep convection parameterization in TJ2016 (see their Fig. 2c). It seems that the modified dynamical core enhances the updrafts and provides enhanced precipitation rate in the tropics, which partially compensates the inefficient vertical transport of the moist air into upper troposphere due to the absence of the deep convection parameterization in present idealized moist Held–Suarez test. In addition, Figs. 17b and 17d highlight the close resemblance between the time-mean zonal-mean specific humidity distributions in both simulations.

Download

Larger image

PowerPoint slide

Fig  17.  As in Fig. 16, but for (a, c) vertical velocity and (b, d) specific humidity. Note that the latitude range for vertical velocity is 15°S–15°N.

In Fig. 18, we use the vertically integrated quantities of full mass $m$, total column water mass ${m_{\rm q}}$, and dry mass ${m_{\rm d}}$ defined by Eqs. (14)–(16) to further explore the improved performance of the modified dynamical core in mass conservation. For this moist Held–Suarez test, total column water mass obviously decreases during the first 5–20 days and then remains basically unchanged. Once again, for the 1200-day long-term integration, one can see that the dry (full) mass increases linearly with time in the modified (original) dynamical core, which is mainly attributed to the non-conservative numerical treatment. However, if we focus on the first 20 days (inset in Fig. 18), one can still see that the full mass is basically conserved and therefore there exists a false anticorrelation between the change of dry mass and total water mass for the original dynamical core (Fig. 18b); while for the modified dynamical core, the dry mass is basically conserved and therefore the change of full mass is now positively correlated with the change of total water mass (Fig. 18a).

Download

Larger image

PowerPoint slide

Fig  18.  Time changes of global mean full mass, total column water mass, and dry mass from day 0 during the 1200-day integration simulated with (a) the modified dynamical core and (b) the original one. Other details are the same as Fig. 10. The inset is an expanded view of the first 20 days, during which the results are not yet contaminated too much by the non-conservative numerical treatment.

5.   Summary

As a comparison to Peng et al. (2019), this paper focuses on the performance of the newly developed nonhydrostatic global spectral dynamical core based on the dry-mass vertical coordinate in applications on the sphere.

First, some 3D dry benchmark tests were conducted to validate the present nonhydrostatic global spectral dynamical core, including steady state, baroclinic waves, mountain waves in non-sheared and sheared background flows, and dry Held–Suarez test. It should be noted that the mountain wave tests are conducted on a reduced-radius sphere, which highlights the nonhydrostatic effects. The steady state test shows that the nonhydrostatic dynamical core can maintain the balanced steady state very well, even for the simulations with very large time step of 3600 s or with very coarse horizontal resolution T_L39, and the ${l_2}$ errors exhibit nearly flat profiles without almost increase over time; the baroclinic wave test shows that the simulated baroclinic wave life cycle is almost identical to that reported in JW2006; the mountain wave results exhibit good quantitative agreement with the corresponding linear analytic solutions and the numerical solutions produced by other models (e.g., MPAS and NH-IFS); and the dry Held–Suarez test shows that the nonhydrostatic dynamical core perfectly reproduces the zonal-mean climatologies from CAM-EUL simulation at a comparable horizontal resolution, although the global-mean mass increases linearly with time, due to the non-conservative numerical treatment (i.e., the semi-Lagrangian advection scheme; Wong et al., 2013). The results from these 3D idealized test cases demonstrate that the present dynamical core is accurate and robust in applications without moisture on the sphere, especially for addressing the nonhydrostatic effects.

Second, three additional moist test cases with simple physics were conducted to further compare the performance of the modified dynamical core and the original Aladin-NH like one. The essential differences between two dynamical cores lies on whether the precipitation mass sink effect is included correctly or not, and thus, the variance between the corresponding simulations depends on the intensity of precipitation. As a case with relatively weak precipitation, the moist baroclinic wave test indicates that the baroclinic instability solutions obtained with these two dynamical cores show no significant differences. However, for testing idealized TC, typically a heavily precipitating system, there are obvious differences between the dynamical cores in both the intensity and location of TC. The modified one tends to obtain simulated results with lower pressure and stronger wind speeds, which are consistent with previous studies (Lackmann and Yablonsky, 2004). Furthermore, the modified one also tends to obtain simulated results with faster northward movement, which we have also found for the splitting supercell test in Peng et al. (2019). Overall, the idealized TC simulated by the modified dynamical core is much closer to the results from the MPAS model, although there is still a slight displacement in the cyclone location.

Considering the uncertainties of idealized TCs simulated by those participating grid models in DCMIP2016, our present results can also be regarded as a supplement to the verifiable database contained in DCMIP2016 for this particular test case. The nonlinear dynamics–physics moisture feedbacks in long-term integration are further explored by the moist Aqua-planet variant of the Held–Suarez test. The comparison results further confirm that the modified dynamical core enhances the updrafts and provides enhanced precipitation rate in the tropics, which partially compensates the inefficient vertical transport of the moist air into upper troposphere due to the absence of the deep convection parameterization in this moist idealized test case. As to the improved performance in mass conservation, if we focus on the first 20-day integration in the moist Held–Suarez test, one can still see that the full mass is basically conserved and therefore there exists a false anticorrelation between the change of dry mass and total water mass for the original dynamical core; while for the modified dynamical core, the dry mass is basically conserved and therefore the change of full mass is now positively correlated with the change of total water mass. Further work should advance the modified dynamical core to full-physics global nonhydrostatic numerical weather prediction application.

In addition, the fluxes of precipitation as well as sedimentation not only transfer mass but also do momentum and energy (e.g., Ooyama, 2001; Ding and Pierrehumbert, 2016; Li and Chen, 2019). All of these transfers are important for simulating condensable-rich planetary atmospheres. However, as the first step to deal with such issues, only the mass transfer by the precipitation flux has been taken into account in present nonhydrostatic dynamical core. In the follow-up work, the momentum and energy transfer associated with the precipitation flux should be further considered through modifying the equations for both the momentum and the enthalpy (temperature) consistently.

Appendix

List of Symbols in Table 1

${\pi _{\rm d}}$, $\pi$—hydrostatic-pressure of dry and moist air;

${\pi _{{\rm ds}}}$, ${\pi _{\rm s}}$—surface hydrostatic-pressure of dry and moist air;

${\dot \pi _{\rm{d}}}$, $\dot \pi$—Lagrangian derivative of dry and moist hydrostatic pressure;

${\eta _{\rm{d}}}$, $\eta$—hybrid terrain-following vertical coordinate based on the mass of dry and moist air;

$A$, $B$—functions used in the definition of the vertical coordinate and are given by predefined values. Here, these predefined values of $A$ and $B$ coefficients for different model levels are directly obtained from https://confluence.ecmwf.int/display/UDOC/Model+level+definitions;

$p$— full pressure (vapor plus dry air);

${\rho _{\rm{d}}}$, $\rho$—density of dry and moist air;

${m_{\rm{d}}}$, $m$—vertical metric factor;

$\hat q$, $dl$—nonhydrostatic prognostic variables for pressure deviation and vertical momentum equations;

$d$, $\chi$—(pseudo) vertical divergence and so-called “X- term”;

${\boldsymbol u}$, $w$—horizontal wind velocity vector and vertical wind velocity component;

${T_{\rm{m}}}$, $T$—modified temperature and original temperature;

${q_j}\left( {j = {\rm{v,c,r}}, \cdots } \right)$—mass mixing ratios of water vapor, cloud water, rainwater, etc.;

$\gamma ={\left(1+{q}_{\rm{v}}+{q}_{\rm{c}}+{q}_{\rm{r}}+\cdots \right)}^{-1}$—ratio of dry air density to moist air density;

${R_{\rm{d}}}$, $R$—gas constant of dry and moist air;

${R_{\rm{v}}}$—gas constant of water vapor;

${c_p}$, ${c_v}$—specific heat of moist air at constant pressure and constant volume;

$Q$—physical contributions to $T$;

${F_\psi }$—physical contributions to $\psi$, such as ${\boldsymbol u}$, $w$, ${q_j}$;

$\varepsilon = {{{R_{\rm{v}}}} \mathord{\left/ {\vphantom {{{R_{\rm{v}}}} {{R_{\rm{d}}}}}} \right. } {{R_{\rm{d}}}}}$—ratio of water-vapor gas constant to dry-air gas constant;

${H_{\rm{m}}} = \left( {1 + \varepsilon {q_v}} \right){Q \mathord{\left/ {\vphantom {Q {{c_p}}}} \right. } {{c_p}}} + \varepsilon T{F_{{q_v}}}$—combined diabatic contri bution;

${D_3}$—true 3D divergence of the wind vector;

$g$—gravitational acceleration;

$\phi$—geopotential;

$\dot \chi$—Lagrangian derivative of $\chi$.