The Non-hydrostatic Icosahedral Atmospheric Model: description and development

Abstract

This article reviews the development of a global non-hydrostatic model, focusing on the pioneering research of the Non-hydrostatic Icosahedral Atmospheric Model (NICAM). Very high resolution global atmospheric circulation simulations with horizontal mesh spacing of approximately O (km) were conducted using recently developed supercomputers. These types of simulations were conducted with a specifically designed atmospheric global model based on a quasi-uniform grid mesh structure and a non-hydrostatic equation system. This review describes the development of each dynamical and physical component of NICAM, the assimilation strategy and its related models, and provides a scientific overview of NICAM studies conducted to date.

1 Review

1.1 Introduction

Diabatic heating due to the release of latent heat in deep convection is the primary heat source in the atmosphere, and it is interacted with the atmospheric general circulation, especially the tropical large-scale overturning circulations such as the Hadley and Walker circulations. Individual deep convective cells are associated with meso-scale circulations that have a horizontal scale of O (10 km), and an upward convective core, along with a horizontal scale of O (km). Until recently, since the horizontal resolution of the global climate models that have been used for future climate change projections has been O (100 km), such models require the use of cumulus parameterizations in order to incorporate the effects of deep convection instead of by explicitly resolving deep convective circulations. However, it is known that cumulus parameterizations significantly affect the results of climate model simulations and that they are the most ambiguous factor used in climate models (Randall et al. [2003]).

To overcome the above-mentioned cumulus parameterization issue, global non-hydrostatic models that utilize a horizontal mesh interval of O (km) for global atmospheric circulation simulations have been developed. Such models explicitly calculate deep convective circulations over the global domain without using cumulus parameterizations. At the grid-resolvable scale, water vapor is saturated into the liquid or ice phase of water in the upward flow field in order to form clouds and is eventually converted to rain and snow through cloud microphysics processes. In global non-hydrostatic models, clouds are spontaneously organized and the multi-scale structures of convective systems are reproduced over the global domain.

The Non-hydrostatic Icosahedral Atmospheric Model (NICAM) (Tomita and Satoh [2004]; Satoh et al. [2008]; Satoh [2013]) was first designed to be run with a horizontal mesh size approximately 3.5 km over the global domain by using the Earth Simulator (http://www.jamstec.go.jp/es/en/) which was launched by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) in 2002. NICAM uses an icosahedral grid, as shown in Figure 1. Higher resolution grids are recursively subdivided from a coarser resolution grid. Hereinafter, we will refer to the grid division level as the g-level. The number of points, arcs, and triangles of the icosahedral grids with g-level l are given as follows:

$$N_P=10n^2+2=10\times 2^{2l}+2\text{,}$$

$$N_A=30n^2=30\times 2^{2l}\text{,}$$

$$N_T=20n^2=20\times 2^{2l}\text{,}$$

Figure 1

Icosahedral grids. Icosahedral grids for grid division levels of 0, 1, and 2 (top, from left to right), and 3, 4, and 5 (bottom, from left to right).

where n = 2^l. We define the average area of the triangles $\overline{A}$ and the average grid interval $\overline{\Delta }$ as follows:

$$\overline{A}=\frac{4\pi R^2}{N_T}=\frac{\pi R^2}{5\times 2^{2l}}\text{,}$$

$$\overline{\Delta }=\sqrt{2\overline{A}}=\sqrt{\frac{2\pi }{5}}\frac{R}{2^l}\text{,}$$

where an Earth radius of R = 6,371.22 km is used. The values for each g-level are listed in Table 1.

Table 1 Number of points, arcs, triangles, and average grid interval of the icosahedral grids

NICAM has been shown to reproduce a realistic multi-scale cloud structure from a meso-scale to a planetary-scale cloud organization that is associated with the Madden-Julian Oscillation (MJO) (Madden and Julian [1971], [1972]) at a g-level between 9 ($\overline{\Delta }$ = 14 km) and 11 ($\overline{\Delta }$ = 3.5 km) (Tomita et al. [2005]; Miura et al. [2007b]). By using the K computer, which is installed at the RIKEN Advanced Institute for Computational Science (AICS) in Kobe, Japan (http://www.nsc.riken.jp/index-eng.html), the resolution of NICAM has been recently increased to the subkilometer level, and it was shown that the deep convective core is more realistically resolved by using a g-level 13 ($\overline{\Delta }$ = 870 m) mesh simulation (Miyamoto et al. [2013]; Figure 2). In the current study, various experiments including decadal continuous experiments and case sweep experiments (Miyakawa et al. [2014]) were also conducted at g-levels between 9 and 11. Future projection studies such as those investigating changes in clouds and tropical cyclone activities are also investigated. The results are interpreted based on the more physically based cloud microphysics processes without the ambiguities of cumulus parameterizations.

Figure 2

Cloud distribution simulated by the NICAM 870 m grid spacing experiment for 6:00 UTC 25 August 2012 (Miyamoto et al. [2013] ).

This article describes the current development status, design, and concepts behind the individual components of NICAM. First, the background and an overview of studies related to the global non-hydrostatic model are reviewed. Then, the history of NICAM development and the scientific outcomes are summarized in the ‘Scientific overview of NICAM’ section, while the NICAM computational design is described in the ‘Design, structure, development, and timeline’ section. The two sections that follow describe the dynamical and physical components of NICAM, and their respective subsections describe each component of the physical processes. Next, the assimilation strategy is described. Finally, various NICAM usages that have been developed by modifying the original NICAM geometry are presented.

1.1.1 Review of global non-hydrostatic models

Thanks to the significant advances in high-performance computers over the last decade, global atmospheric simulations with a horizontal resolution of O (10 km) can be achieved (Ohfuchi et al. [2004]; Mizuta et al. [2005]; Kinter et al. [2013]; Wedi [2014]). At enhanced horizontal resolution of less than 10 km, traditional atmospheric general circulation models (AGCMs) encounter fundamental difficulties in their dynamic framework formulation as well as in their computational efficiency. As their resolutions increase, AGCMs can capture flow features with comparable scales of motion in the horizontal and vertical directions (such as deep convection and fine-scale gravity waves) that can invalidate the hydrostatic approximation. Although deep convective systems in the tropics play key roles in global atmospheric circulations, they have not been directly resolved by AGCMs, and their effects have only been considered in a parameterized form (Arakawa [2004]). The effects of fine-scale gravity waves that are parameterized as gravity wave drag in AGCMs can be captured as the resolution increases, but the propagation of such gravity waves will be incorrectly calculated unless the non-hydrostatic effect is taken into consideration (Iwasaki et al. [1989]).

Regarding the numerical algorithm, many of the existing AGCMs employ the spectral transform method to represent spherical fields, and it has been pointed out that spectral transforms become increasingly inefficient for high-performance computing as the horizontal resolution increases (e.g., Stuhne and Peltier [1996]; Taylor et al. [1997]; Randall et al. [2000]; Satoh et al. [2005]; Cheong [2006]; Tomita et al. [2008]; Wedi [2014]). Another problem that occurs during computation on a massively parallel computer is that the spectral transform method requires extensive data movement between computer nodes. Although the double Fourier transformation method has been proposed as an alternative (e.g., Cheong [2006]), this method still requires global communication between computer nodes.

To increase the horizontal resolution beyond O (10 km) in a global atmospheric model, the governing equations and numerical algorithms must be reconsidered. More specifically, the governing equations must be non-hydrostatic, and the grid point method replaces the spectral method. As for the familiar latitude-longitude grid (lat-lon grid), however, the grid spacing near the poles becomes drastically reduced as the horizontal resolution is increased, which means that reduced grids are generally required to avoid severe time interval restrictions for the Courant-Friedrichs-Lewy (CFL) condition. In principle, the semi-Lagrangian, semi-implicit (SLSI) approach (cf. Laprise [2008]; Staniforth and Wood [2008]) could be employed to overcome the requirement of the CFL condition. Several authors (Semazzi et al. [1995]; Cullen et al. [1997]; Qian et al. [1998]; Côté et al. [1998]; Yeh et al. [2002]; Davies et al. [2005]; Wedi and Smolarkiewicz [2009]; Wood et al. [2013]) have used the lat-lon grid to solve a set of non-hydrostatic equations using the SLSI approach in order to acquire a larger time interval for integration. However, it is unclear how effective the elliptic solvers, which were developed for SLSI schemes, would be for ultra-high-resolution calculations, even though it has been recently proven that a multi-grid approach is an ideal solver for massively parallel computers (Heikes et al. [2013]).

The pole problem can be overcome by using grid systems with quasi-homogeneous grids over the sphere. One such grid is the icosahedral grid, which is currently one of the major grid systems used for high-resolution global atmospheric modeling. Primitive (hydrostatic) equation global models using icosahedral grids have been developed at Colorado State University for climate modeling (CSU AGCM) (Ringler et al. [2000]), Deutscher Wetterdienst for numerical prediction modeling (GME) (Majewski et al. [2002]), the National Oceanic and Atmospheric Administration (NOAA) for the Flow-following finite-volume Icosahedral Model (FIM) (http://fim.noaa.gov/), and the Laboratoire de Météorologie Dynamique (DYNAMICO) (http://www.lmd.polytechnique.fr/~dubos/DYNAMICO/). Global non-hydrostatic models using icosahedral grids are also being developed by several international groups including the geodesic grid model at Colorado State University (UZIM; http://www.cmmap.org/research/models.html), the Icosahedral Non-hydrostatic model (ICON) at the Deutscher Wetterdienst and the Max Planck Institute for Meteorology (http://www.mpimet.mpg.de/en/science/models/icon.html; Zängl et al, [2014]), the Model for Prediction Across Scales (MPAS) at the National Center for Atmospheric Research (Skamarock et al. [2012]), and the Non-hydrostatic Icosahedral Model (NIM) at NOAA (http://www.esrl.noaa.gov/gsd/ab/ac/GPU_Parallelization_NIM.html). Entries for all these models can also be found at https://www.earthsystemcog.org/projects/dcmip-2012/. As for other types of grid models, cubic grids are also candidates for high-resolution atmospheric models (McGregor [1996]; Lin [2004]; Putman and Suarez [2011]).

In Japan, the icosahedral atmospheric model has been created by using a non-hydrostatic system, i.e., NICAM (Tomita and Satoh [2004]; Satoh et al. [2008]; Satoh [2013]; http://nicam.jp/). Development of NICAM began around 2000. Since high-resolution modeling has now entered the mainstream worldwide, NICAM has joined a number of international high-resolution numerical modeling projects, such as the Athena Project (Kinter et al. [2013]) with the Integrated Forecast System (IFS) (Jung et al. [2012]) and the Icosahedral-grid Models for Exascale Earth system simulations (ICOMEX) (Zängl et al. [2011]) with ICON, MPAS, and DYNAMICO.

1.1.2 Scientific overview of NICAM

While more than a decade has passed since the development of NICAM began, the first milestone was reached in 2004 when Tomita and Satoh ([2004]) began finalizing the dynamical core of the model. Two unique characteristics of the NICAM dynamical core are the use of the spring dynamics to construct a modified icosahedral grid (Tomita et al. [2001], [2002]) and the use of a flux form non-hydrostatic system which guarantees conservation of mass and energy (Satoh [2002], [2003]). Tomita et al. ([2001], [2002]) constructed NICAM's icosahedral grids by improving the numerical accuracy of differential operators and using spring dynamics to improve the homogeneity of the grid system. An improved version of the spring dynamics method, which is more homogeneous and applicable even for higher resolution, is also now available (Iga and Tomita [2014]). NICAM uses a fully compressible (elastic) non-hydrostatic system to obtain statistically equilibrium states by performing long-timescale simulations. For this purpose, we devised a non-hydrostatic numerical scheme that guarantees the conservation of mass and energy (Satoh [2002], [2003]). Tomita and Satoh ([2004]) implemented this non-hydrostatic scheme in their global model using an icosahedral grid configuration, with which they developed a dynamical core for NICAM. Since the split-explicit time integration scheme (Klemp and Wilhemson [1978]) is used for the horizontal propagation of fast waves and for the implicit treatment of the vertical propagation of fast waves, multi-dimensional elliptic solvers are not required. To extend the original non-hydrostatic scheme to the global domain, the set of equations provided in Satoh ([2002], [2003]) is reformulated for spherical geometry and modified in order to make the equations suitable for an icosahedral grid configuration. The finite volume method is used for numerical discretization, so that the total mass and energy over the domain are conserved. The resulting model is suitable for long-term climate simulations.

NICAM is used for ‘cloud-resolving simulations’ by explicitly resolving convective circulation. This is accomplished by drastically increasing the horizontal resolution using high-power computing systems, such as Earth Simulator and the K computer. After finalizing the development of the dynamical core (Tomita and Satoh [2004]), the first 3.5 km mesh global NICAM simulation was performed on Earth Simulator by Tomita et al. ([2005]) using the aqua planet configuration (Neale and Hoskins [2001]). The results of this simulation clearly show multi-scale convective systems propagating near the equator (Satoh et al. [2005]), which is similar to the propagation of the observed cloud clusters and super cloud cluster structure (Hayashi and Sumi [1986]; Nakazawa [1988]; Takayabu et al. [1999]).

Tomita et al. ([2005]) conducted aqua planet experiments at three horizontal resolutions, g-levels 9, 10, and 11 ($\overline{\Delta }$ = 14, 7, 3.5 km), using the same physical schemes without cumulus parameterization. In their results, the multi-scale convective structure along the equatorial zone and the diurnal cycle of convective precipitation were reproduced similarly in all simulations, although the propagation speed along the equator and the phase lag of the diurnal cycle depends on the resolution. This result motivated us to use a combination of different horizontal resolution experiments with the same physical schemes when studying convective properties simulated by NICAM.

The results of the NICAM aqua planet experiment were analyzed intensively by Nasuno et al. ([2007], [2008]), Nasuno ([2008]), and Nasuno and Satoh ([2011a], [b]), which showed the roles of the equatorial convective systems embedded in the multi-scale structure. Mapes et al. ([2008]) also investigated the persistent structure in the tropics in terms of the predictability of disturbances. The diurnal cycle of tropical convective precipitation was clearly simulated with the dependency of the phase lag on the horizontal resolution (Tomita et al. [2005]). The semi-diurnal cycle was also captured. This was further investigated by Yasunaga et al. ([2013]) by using the NICAM aqua planet experiment. The detailed structures of convective systems near the tropopause were analyzed by Kubokawa et al. ([2010]).

The aqua planet configuration is used to test cloud changes that are related to increased surface temperatures in order to imitate global warming conditions. Miura et al. ([2005]) compared the cloud cover responses between NICAM and the Model for Interdisciplinary Research on Climate (MIROC) (Hasumi and Emori [2004]), which is an ordinary resolution climate model (denoted as CCSR/NIES/ FRCGC in Miura et al. [2005]), under the aqua planet condition. The comparison results show that the cloud cover simulated by NICAM increases in high-latitude regions, while that by MIROC decreases. Such a contrasting response to cloud cover is very interesting and should be analyzed in more detail, as will be described later in this section. The change in the meridional distribution of relative humidity is shown and compared to the other model results by Sherwood et al. ([2010]).

The aqua planet experiments of NICAM were compared to the other model results in Blackburn et al. ([2013]), and additional aqua planet experiments were conducted by Yoshizaki et al. ([2012a], [b]) in order to investigate the multiple-scale convective structure along the equator.

A NICAM simulation that includes a realistic configuration with a land and sea contrast was performed by Miura et al. ([2007a]). The results of this simulation show not only realistically simulated tropical cyclogenesis but also a clear dependency of convective organization on a planetary boundary layer (PBL) scheme, specifically, water vapor transport. Using the simulated dataset, the diurnal variation of the convective systems over the Tibetan Plateau was analyzed by Sato et al. ([2007], [2008]).

A successful simulation of a realistic MJO event was presented by Miura et al. ([2007b]), in which a 1 month integration with a 7 km grid spacing and a 1 week integration with a 3.5 km grid spacing were performed. A large-scale cloud organization, its eastward propagation, and the multiple-scale structure of convective systems were realistically reproduced (Nasuno et al. [2009]). The similarity of the simulation to the observed MJO was further analyzed by Liu et al. ([2009]). In addition, multiple tropical cyclones were realistically generated from the active cloud areas of the MJO. Particularly, Fudeyasu et al. ([2008], [2010a], [2010b]) analyzed the evolution of a tropical storm in detail, which was generated from the MJO 2 weeks after the initial condition.

Thanks to the high-resolution numerical simulations performed by Miura et al. ([2007b]), simulation results can be analyzed or used in various ways. The dataset consists of 7 km mesh and 14 km mesh grid data for 1 month from 15 December 2006 and 3.5 km mesh grid data for 1 week from 25 December 2006. The diurnal cycle of tropical convective systems was intensively analyzed by Sato et al. ([2009]). Convective momentum transports of the MJO were analyzed by Miyakawa et al. ([2012]). The dataset was compared with the satellite observations of the Tropical Rainfall Measuring Mission (TRMM) and CloudSat (Masunaga et al. [2008]), the geostationary satellite, the Multi-functional Transport Satellite (MTSAT-1R) (Inoue et al. [2008]), and CloudSat and Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) (Inoue et al. [2010]; Satoh et al. [2010]; Ham et al. [2013]), and with the estimation of cloud radiative forcing (Sohn et al. [2010]). The dataset was also used to develop a cloud parameterization for coarser resolution AGCMs (Watanabe et al. [2009]). The convective structure near the tropopause was also analyzed by Kubokawa et al. ([2012]). The energy spectrum and similar frequency spectra of the high-resolution data were analyzed by Terasaki et al. ([2009]) and Tsuchiya et al. ([2011]).

Following Tomita et al. ([2005]) and Miura et al. ([2007b]), studies of intra-seasonal variability (ISV) including MJOs and boreal summer intra-seasonal oscillations (BSISO) (Fu and Wang [2004]) are unique areas of NICAM studies. Various aspects of the MJO and tropical disturbances simulated by Tomita et al. ([2005]) and Miura et al. ([2007b]) are being analyzed as mentioned above. NICAM studies are also collaborated with observational field campaigns of MJOs. Miura et al. ([2009]) examined an onset of MJO observed at the Mirai Indian Ocean cruise for the Study of the MJO-convection Onset (MISMO) field campaign (Yoneyama et al. [2008]), suggesting a role of the sea surface temperature (SST) gradient in the western Indian ocean. Also, for the Cooperative Indian Ocean Experiment on Intraseasonal Variability in the Year 2011/Dynamics of the MJO (CINDY2011/DYNAMO) campaign (Yoneyama et al. [2013]), Nasuno ([2013]) summarized a result from quasi-real-time weather forecasting experiments that showed simulated characteristics of MJOs during the CINDY2011/DYNAMO period in the boreal winter between 2011 and 2012. The quasi-real forecasting system using the stretch NICAM (Tomita [2008a]; see ‘Models related to NICAM’ section) was developed particularly for use in collaboration with field experiments, as introduced by Oouchi et al. ([2012]). More statistical results for MJO simulations are shown by Miyakawa et al. ([2014]), who reported a good performance of MJO predictability by the hind cast approach with NICAM. BSISO simulations were analyzed by Oouchi et al. ([2009a], [2014]), Taniguchi et al. ([2010]), Yanase et al. ([2010b]), and Satoh et al. ([2012a]). In particular, they focused on relations between ISV and tropical cyclones, as described below, since ISV modulates large-scale vorticity distributions that generally lead to preferable conditions for cyclogenesis.

Tropical cyclones can also be adequately researched by using NICAM. Tropical cyclogensis and its evolution from the MJO were analyzed for boreal winter cases by Fudeyasu et al. ([2008], [2010a], [2010b]). For boreal summer, tropical cyclogenesis was found to be frequently related to ISV activities (Kikuchi and Wang [2010]). The cyclone Nargis, which caused severe damage to Myanmar in May 2008, was successfully simulated and its relation to the northward propagation of BSISO in the Bay of Bengal was analyzed by Taniguchi et al. ([2010]) and Yanase et al. ([2010a], [2012a]). Oouchi et al. ([2009a]) simulated tropical cyclogenesis in the northwestern Pacific and showed its relation to the MJO. The effects of an equatorial Kelvin wave on an abrupt development of a tropical cyclone were clearly shown by Yanase et al. ([2010b]). Yanase et al. ([2012b]) also analyzed some statistical behaviors of tropical cyclones simulated by NICAM especially over the eastern Pacific. The impacts of a tropical cyclone on the Baiu rainfall were analyzed by Yamaura et al. ([2013]). The preferable conditions of tropical cyclogenesis were analyzed over the Atlantic domain (Satoh et al. [2013]).

Tropical cyclone climatology and possible changes in the future projection of tropical cyclone activities are actively being researched using NICAM. For example, Satoh et al. ([2012a]) analyzed 8 years of boreal summer experiments that were conducted in the Athena Project (Kinter et al. [2013]) in order to show the simulated climatology of tropical cyclones and produced results that show good relations between tropical cyclone activities and MJO activities. Future changes in the tropical cyclone structure were analyzed by Yamada et al. ([2010], [2012]), Yamada and Satoh ([2013]), and Emanuel et al. ([2010]).

NICAM simulations also contribute to the studies on Asian monsoon and tropical convective systems. For example, Oouchi et al. ([2009b]) provided a realistic distribution of precipitation during a boreal summer and showed its variation for diurnal and intra-seasonal timescales. Additionally, the resolution dependency of the diurnal cycle of convective systems was analyzed by Sato et al. ([2009]), and their results were found to be generally consistent with those of the aqua planet experiment (Tomita et al. [2005]); as resolution becomes coarser, the timing of precipitation peak becomes increasingly delayed. Sato et al. ([2009]) also analyzed the diurnal evolution of tropical convective cold pools, while Fujita et al. ([2011]) analyzed the behavior of the diurnal convection over maritime continental regions. Dirmeyer et al. ([2012]) and Noda et al. ([2012]) also analyzed the diurnal cycle of precipitation over the global domain. A mechanism for the diurnal variation of summer precipitation over southern China was investigated by Satoh and Kitao ([2013]) using the stretched version of NICAM (see ‘Models related to NICAM’ section).

The use of a cloud microphysics scheme in NICAM experiments without cumulus parameterization leads to realistic behavior of cloud distribution in the horizontal and vertical directions. Furthermore, as seen in the paragraph below, it might even produce a different behavior of cloud feedback from that obtained by conventional climate models that use cumulus parameterization. Iga et al. ([2007a]) analyzed the climatological distribution of cloud cover and compared it with the data of the International Satellite Cloud Climatology Project (ISCCP) (Rossow and Schiffe [1999]), while in a later study, Iga et al. ([2011]) analyzed upper tropospheric ice clouds and their sensitivity using various model parameters. Noda et al. ([2010]) investigated the importance of the subgrid turbulent process on cloud distributions. Kodama et al. ([2012]) showed the sensitivities of upper clouds to the choice of various cloud microphysics parameters. Satoh and Matsuda ([2009]) investigated how cloud microphysics affects upper clouds under an idealized radiative convective equilibrium condition.

As for the future simulations of cloud change by NICAM, Miura et al. ([2005]) showed the model's different behavior when compared to a coarse resolution AGCM (MIROC). In addition, Collins and Satoh ([2009]) showed that as a response that was different from the existing models, the upper cloud cover might increase, whereas the ice water path decreases under future warming conditions. Satoh et al. ([2012b]) analyzed their results and speculated that a future reduction in convective circulation could lead to a reduction in the ice water path under the warming condition. Tsushima et al. ([2014]) further investigated longwave cloud feedback and showed the possibility of a stronger positive feedback than other coarse resolution AGCMs. Yamada and Satoh ([2013]) analyzed the relation between cloud changes over the global domain and tropical cyclone changes. Kodama et al. ([2014a]) also analyzed the cloud changes associated with storm tracks simulated by the NICAM aqua planet experiments.

Recently, use of the K computer has allowed the horizontal resolution of NICAM to be increased. Miyamoto et al. ([2013]) performed (for the first time) a subkilometer mesh global domain simulation with grid spacing of 870 m (g-level 13) and 96 vertical layers (Figure 2). They showed that at a resolution of less than 1.7 km, convective cores could be resolved by multiple grid cells. The ability to use the K computer introduces a new area of study for NICAM. Numerous MJO case studies are currently being conducted in order to evaluate their statistical predictability (Miyakawa et al. [2014]). More than 10 years of continuous simulations have and are being conducted under present and future conditions in order to investigate future changes to tropical cyclones and other extreme events. The subsequent sections describe the improvements in the simulated climatology that have been made possible by improving the physical processes.

As described above, NICAM is usually run without using cumulus parameterization. In general, the horizontal scale of meso-scale phenomena ranges from a few to several hundred kilometers. It is not readily understood why a model with horizontal mesh spacing of O (10 km) reproduces such meso-scale phenomena without subgrid convective parameterization. It is speculated that for relatively large-scale meso-scale convective systems, the statistical effects of circulation resolved by grids behave similarly to nature. Figure 3 shows resolution dependency of upper cloud distributions in the range of $\overline{\Delta }$ = 28, 14, 7, 3.5, 1.7, 0.87 km for the NICAM simulations conducted by Miyamoto et al. ([2013]). Although in that study, it was analyzed that the convective core becomes resolved at $\overline{\Delta }$ = 1.7 km. Figure 3 shows that the large-scale structure of upper clouds is nearly unchanged for these resolutions. This result implies that the usefulness of a model without cumulus parameterization depends on research targets of meteorological phenomena. Large-scale organized convective structures such as MJOs are, in general, well reproduced without using cumulus parameterization, even at the coarser horizontal resolutions around O (10 km) as also indicated by Holloway et al. ([2013]), while the precise simulation of details of meso-scale convective structure needs much finer resolutions (Bryan et al. [2003]).

Figure 3

Resolution dependency of cloud structure simulated by NICAM between 28 km and 870 m grid spacing experiments. Cloud distributions simulated by NICAM between 28 km and 870 m grid spacing experiments for 6:00 UTC 25 August 2012, shown by the outgoing longwave radiation at the top of the atmosphere. The unit is W m⁻².

1.1.3 Design, structure, development, and timeline

A timeline of the development of NICAM is summarized in Table 2. As was previously mentioned, the development of NICAM started in 2000 at JAMSTEC. Through the construction of icosahedral grid systems (Tomita et al. [2001], [2002]), a shallow water model (SWM) (Tomita et al. [2001], [2002]) was initially tested. At the same time, new non-hydrostatic numerical schemes were being developed (Satoh [2002], [2003]). By combining the two approaches, a three-dimensional (3D) dynamical core was eventually developed (Tomita and Satoh [2004]). The resulting 3.5 km mesh high-resolution dynamical core experiments were examined by Iga et al. ([2007b]). Then, full-physics experiments have been underway since 2004. Furthermore, the dynamical schemes, computational stability, and a tracer advection scheme (Miura [2007]; Niwa et al. [2011a]) have improved. For physical processes, we have actively introduced more sophisticated schemes, especially for cloud microphysics (Grabowski [1998]; Tomita [2008b]; Seiki and Nakajima [2014]) and PBL schemes (Nakanishi and Niino [2006], [2009]; Noda et al. [2010]). An accurate semi-Lagrangian scheme for sedimentation processes was developed for cloud microphysics schemes (Xiao et al. [2003]). Aerosol processes (Suzuki et al. [2008]) and atmospheric chemistry processes have also been introduced. Each component of the abovementioned physical schemes will be described in later sections of this review. Several NICAM-based systems are also currently under development, including atmosphere-ocean coupling and ensemble-based data assimilation systems, and downscaling methods for regional scale experiments have been developed for more advanced uses.

Table 2 Timeline of the development of NICAM

Since high-resolution simulations require large computational resources, the efficiency of state-of-the art supercomputers needs to be maximized. In general, weather/climate models use a large number of physical variables during computation, so the size of 3D arrays significantly increases with increases in spatial resolution. Large-sized, various, frequent data output is also needed for simulations. Taken together, these factors make it clear that weather/climate models require faster data throughput in all layers of a computer system. NICAM has been tactically designed to achieve efficient throughput at the 1) memory layer, 2) network layer, and 3) file input/output (I/O) layer.

Since numerous simple arithmetic operations are used for large arrays in NICAM algorithms, memory throughput is the most important aspect governing total computational efficiency. For example, when using Earth Simulator, the ratio of memory bandwidth to the floating-point performance (Byte/Flops ratio; hereinafter B/F ratio) is 4, which makes it easy to execute NICAM with a high degree of efficiency. However, in the past decade, the B/F ratio has decreased during the recent evolution of supercomputers. As a result, when NICAM was implemented on the K computer at B/F = 0.5, numerous optimizations had to be applied to improve the efficiency of computer performance. Saving memory transfers is one effective optimization tactic. To accomplish this, we suppressed the memory copies, reduced the intermediate arrays, and avoided unnecessary zero-filling. These changes had to be performed manually and required labor-intensive efforts. However, as a result of these optimizations, we succeeded in reducing the model execution time by 30%.

To improve network throughput, the finite volume method with an explicit scheme was adopted to minimize the global communications, as described in the section ‘Dynamics’. When 640 nodes on the K computer are used, the network communication of NICAM accounts for about 10 % of the total elapsed time at the case of g-level 9 ($\overline{\Delta }$ = 14 km). This percentage increases as the number of nodes increases. Kodama et al. ([2014b]) developed an algorithm for optimally assigning of nodes in a network topology such as a torus/mesh.

For the file I/O throughput, NICAM adopts a distributed file I/O. If the output data are aggregated into a single file at the time of output, it would be easy to handle the file. However, NICAM cannot maximize the I/O throughput using such a single output file. The data of each region on the icosahedral grid of NICAM are output to a file at the node, and the icosahedral grid is converted into a lat-lon grid during the post-process. Generally, we use a simple bilinear interpolation for the post-process, although we have an option of another interpolation to keep conservation of area integration. A recently developed post-process program can run in parallel with the main model using a coupler (see ‘Coupler’ section). However, an increase in file size due to the higher resolution is a current issue, and we are planning to change the file type to a compressible standard format, such as Network Common Data Form version 4 (NetCDF4) (http://www.unidata.ucar.edu/software/netcdf/). Recently, a benchmark run on the K computer has achieved 10% of the peak performance using five nodes (40 cores), and usage has been verified with 81,920 nodes (655,360 cores). Furthermore, a subkilometer mesh global domain simulation (Miyamoto et al. [2013]) was executed with a performance of 230 Tflops in double precision for 68 billion grid cells while using 20,480 nodes of the K computer.

Further information on the NICAM development team and the initial developments can be found at http://nicam.jp. NICAM team members consist of researchers belonging to several research institutions in Japan. The source code is primarily written in Fortran 95 and follows the Meteorological Research Institute/Japan Meteorological Agency (MRI/JMA) standard coding rule (Muroi et al. [2002]). The program is parallelized using Message Passing Interface (MPI) process parallelization via automatic compiler thread parallelization. OpenMP is not actively used as the current program. Each component is managed as one or more modules. In order to maintain readability of the source code, C preprocessor macros are not used. The source code has been managed at repositories using the Concurrent Versions System (CVS), and the NICAM development team is currently moving to multi-branch development using the Git distributed revision control and source code management system. Due to its flexible portability, NICAM can run on various computer platforms, such as Earth Simulator, the K computer, IBM AIX, Linux, and Mac OSX. We have also tested several combinations of MPI libraries and Fortran compilers.

NICAM is also used as a benchmark program or a test platform for new architectures. For this purpose, we have released the dynamical core of NICAM as open source software (NICAM-DC), which is distributed under the terms of the Berkeley Software Development (BSD) 2-Clause License and is available for download (http://scale.aics.riken.jp/nicamdc/). The NICAM development team also participates in a working group for the common base library environment of weather and climate models and in cooperative development efforts with other modeling groups in Japan. In the field of computational science, NICAM has made significant contributions as a benchmark program and is one of the main application programs used in Exascale computing research.

1.2 Dynamics

1.2.1 Grid configuration and advection scheme

As reviewed in the ‘Introduction’ section, the use of grid models with finite-difference and finite-volume methods for weather and climate models was reconsidered as massively parallel computer architectures entered widespread use. An icosahedral grid is a globally quasi-uniform and is possibly the most suitable base for the development of an atmospheric model that can be run on a supercomputer. Additionally, the adoption of an icosahedral grid for the horizontal grid of NICAM was partially motivated by the results of Heikes and Randall ([1995]), in which a two-dimensional (2D) shallow water system model was developed using the Z-grid arrangement of the prognostic variables (Randall [1994]; Figure 4d), which followed Masuda and Ohnishi ([1986]). The Z-grid model, however, had a disadvantage in that the model needed to diagnose the velocity fields from the vorticity and divergence fields by solving two Poisson equations in each time step. Furthermore, it was not obvious that an efficient Poisson solver would become available for use on the massively parallel computers, at least when the NICAM development began.

Figure 4

Arrangements of variables on the hexagonal grid. Arrangements of prognostic variables on the hexagonal grid for the (a) A-grid, (b) B/ZM-grid, (c) C-grid, and (d) Z-grid staggering. The symbols ρ, T, u, v, v _⊥, δ, and ζ refer to density, temperature, zonal components, and meridional components of the flow velocity, respectively, that is defined (a) at the center or (b) at the vertex of the hexagon, normal component of the flow velocity on each edge, divergence of the horizontal flow, and vertical component of vorticity.

Tomita et al. ([2001]) did not follow any previous study (such as Heikes and Randall [1995]; Masuda and Ohnishi [1986]) and instead examined a new method by using the Arakawa A-grid (Arakawa and Lamb [1977]), which located prognostic variables (fluid depth and velocity) at the same nodes (Figure 4a) and thus did not require Poisson solvers. The accuracy and stability of the shallow water model were improved after applying a grid optimization that is now commonly known as ‘spring dynamics’ (Tomita et al. [2002]), and the physical performance of the model was evaluated using the standard tests provided in Williamson et al. ([1992]). The methods used to discretize the horizontal gradient, divergence, rotation, and Laplacian operators in NICAM are similar to those described in Tomita et al. ([2001]) and Tomita and Satoh ([2004]), except for some minor updates in the passive tracer transport (Miura [2007]; Niwa et al. [2011a]) and damping operators.

When the model was first developed by Tomita et al. ([2005]), NICAM was capable of producing a global simulation with horizontal grid spacing of 3.5 km. It was also the first model to produce a global simulation with horizontal grid spacing of less than 1 km (Miyamoto et al. [2013]). These characteristics indicate that the A-grid arrangement is a good choice due to its simplicity and the relative ease with which it attains high computational performance on Earth Simulator and the K computer. It should be noted that Tomita et al. ([2001]) provided the first example of the shallow water model on an icosahedral grid, which uses the finite-volume method and adopts velocity components as prognostic variables. However, it is well known that the A-grid model is unsuitable for simulations of inertia-gravity waves and geostrophic adjustments and that the model suffers from grid-scale noise (Randall [1994]). Therefore, other options should be considered for the horizontal variable arrangement, such as those described below.

The disadvantages of the B/ZM-grid (Figure 4b) and C-grid (Figure 4c) arrangements, which were pointed out by Ničkovic et al. ([2002]), on the hexagonal grid were partially solved by Ringler and Randall ([2002]) and Thuburn et al. ([2009]), respectively. Although problems with the computational modes remain due to the mismatch of the degrees of freedom between the fluid depth/pressure and horizontal velocity, the B/ZM-grid and C-grid not only allow a more realistic representation of inertia-gravity waves and geostrophic adjustments than the A-grid but also are free from the checkerboard pattern. Currently, we are investigating the B/ZM-grid as a candidate to replace the A-grid, and the difficulty in managing the computational mode has been largely eliminated by revising the operators of Ringler and Randall ([2002]). This activity is ongoing, and details of the improvements will be reported in forthcoming papers. The update of the dynamical core will also include the newly developed transport algorithms of Miura and Skamarock ([2013]) and Miura ([2013]).

1.2.2 Vertical resolution issues

NICAM adopts the Lorenz grid and z* (terrain-following) coordinates (Gal-Chen and Somerville [1975]) as a vertical coordinate system (Figure 3 of Satoh et al. [2008]). Density, horizontal velocity, and internal energy are defined at the full levels, and vertical velocity is defined at the half levels. In NICAM, z* is generally specified so as to linearly depend on the geometric height z between the model surface and top (the linear-z* coordinate). Using the standard vertical levels Z _k , geometric height z _X,k is formulated as

$$z_{X,k}=Z_k+\left(1-\frac{Z_k}{H_{\mathrm{top}}}\right)Z_{X,\mathit{sfc}}\text{,}$$

where k is an index for the vertical level, X represents the location of the horizontal grid, and H _top and Z _X,sfc are the heights of the model top and bottom, respectively. In general, H _top = 40 km is used with 40 vertical levels. For standard cases, the lowermost full level is located at 81 m. The vertical interval of the half-level ΔZ _k increases with height, from ΔZ _k = 170 m at the bottom, ΔZ _k  < 1 km below 10 km, and ΔZ _k  ~ 3 km near the model top.

In some situations, the hybrid-z* coordinate (Simmons and Burridge [1981]) is employed to reduce the amount of pressure gradient force numerical errors that arise from topography. It is formulated as

$$z_{X,k}={\mathit{Z}}_k+\left(\frac{sinh\frac{{\mathrm{H}}_{\mathrm{top}}-Z_k}{H_{\mathrm{efold}}}}{sinh\frac{H_{\mathrm{top}}}{{\mathrm{H}}_{\mathrm{efold}}}}\right)Z_{X,\mathit{sfc}}\text{,}$$

where H _{e fold} is a height scale parameter (e.g., H _{e fold} = 10 km). In the hybrid-z* coordinate, intervals of the vertical level approach those of the standard vertical level for Z _k ≫ H _{e fold}, whereas in the linear-z* coordinate, it is constant with height even near the model top. This hybrid-z* coordinate advantage is obtained at the expense of packing more vertical levels just above the higher mountain, which may cause the vertical CFL condition to become a serious issue for the hybrid-z* coordinate. Klemp ([2011]) also proposed a hybrid approach, in which the coordinate surfaces are smoothed with height in order to remove smaller scale terrain structures from the surfaces.

Recently, we constructed a new configuration for the vertical levels of NICAM in order to perform simulations with finer vertical resolutions in the troposphere, stratosphere, and mesosphere levels. The development of this ‘middle atmosphere NICAM’ is motivated by the results of recent climate studies in which the tropospheric circulations were found to be affected by the stratosphere. It should be noted that the inclusion of the middle atmosphere into a climate model is a worldwide trend in the climate model community, since it improves the reproducibility of the large-scale tropospheric circulation (e.g., Shindell et al. [1999]). Currently, the middle atmosphere NICAM employs the hybrid-z* coordinate, and its model top height is set to 80 km. The vertical intervals in the lower atmosphere are the same as the standard vertical coordinate. As the vertical interval increases with height and reaches a specific criteria (Δz _max), the vertical level is set with the uniform interval of Δz _max (except near the model top) in order to appropriately simulate the propagation of atmospheric gravity waves. We have tested Δz _max = 2, 1, and 0.5 km with a horizontal resolution of g-level 7 ($\overline{\Delta }$ = 120 km), as shown in Figure 5. As the vertical resolution increases, the reproducibility of the zonal mean temperature and zonal wind tends to improve in the upper troposphere and middle atmosphere. Interestingly, westerly winds around the equatorial lower troposphere in June 2004 area only well simulated when the vertical resolution is set to 1 km or less. NICAM has an option that allows the deep atmosphere configuration (Tomita and Satoh [2004]) to be switched. For the middle atmosphere simulations, the choice of the deep atmosphere is considered appropriate.

Figure 5

Dependency of the simulated NICAM climatology on the vertical resolution. Dependency on the vertical resolution of the (upper) zonal mean temperature (K), and (lower) zonal mean zonal wind (m s^-1). From left to right: JRA25 reanalysis (Onogi et al. [2007]), CIRA86 (http://badc.nerc.ac.uk/view/badc.nerc.ac.uk__ATOM__dataent_CIRA), NICAM (Δz_max = 2 km), NICAM (Δz_max = 1 km), and NICAM (Δz_max = 0.5 km). For CIRA86, the climatology for June is shown. All the other panels are the averages for June 2004. Shading denotes deviations from the CIRA86 climatology.

Determining the proper vertical coordinates is an ongoing issue for high-resolution non-hydrostatic models. For global non-hydrostatic models in particular, the numerical instability of the terrain-following coordinates becomes more severe because the steep topography over the Andes or the Antarctica is within the simulation domain. For the present simulations, relatively coarser topographic data were used by applying a smoother for the NICAM simulations. For example, for the 3.5 km mesh simulations, 14 km mesh coarsened topographic data were generally used. In addition, a minimum threshold value of the Brunt-Väisälä frequency was specified for the calculations in the PBL schemes, that is and the Mellor-Yamada Nakanishi-Niino (MYNN) scheme, in order to suppress an unusual increase in vertical shears near the surface.

To overcome the problem of terrain-following coordinates, different candidates have been proposed. For example, the possibility of a vertical coordinate that is based on geopotential height, i.e., the z-coordinate, has been investigated. Although this approach is common for ocean models (Adcroft et al. [1997]), it must be understood that the implementation of the z-coordinate to non-hydrostatic atmospheric models is not straightforward (Gallus and Klemp [2000]). The immersed boundary method is well known for computational fluid dynamics (Mittal and Iaccarino [2003]). Various approaches to overcome the problems of the z-coordinate were taken into consideration, such as the cut-cell method or the thin-wall approximation (e.g., Bonaventura [2000]; Steppeler et al. [2002], [2006]; Yamazaki and Satomura [2008], [2010]), have been taken into consideration. We plan to introduce such an option for the z-coordinate of NICAM in the near future.

1.3 Physics

1.3.1 Cloud microphysics schemes

Cloud microphysics can be used to solve the growth equations of a single hydrometeor particle. The related processes are mainly categorized into four components: nucleation processes, phase change, collisional processes, and gravitational sedimentation. These processes have not been calculated explicitly and have been parameterized in conventional AGCMs since the characteristic timescale of these processes and spatial variability of cloud ensembles cannot be resolved using the time-step and grid resolution of GCMs. Global non-hydrostatic models operate by calculating these processes via cloud microphysics schemes in order to explicitly evaluate time evolution of cloud growth over the global domain. The implemented schemes are generally developed for regional non-hydrostatic models. This section introduces cloud microphysics models that are applicable to global non-hydrostatic models. For a more fundamental theoretical basis of cloud microphysics processes, refer to Pruppacher and Klett ([1997]).

Multi-moment bulk method

Hydrometeor particles typically have a radius r from 10 μm to 1 mm and populations from 10² to 10⁹ m⁻³. To manage numerous particles efficiently (in terms of computational cost), cloud microphysics modeling generally assumes that particles with the same radius are homogeneously distributed within a grid cell of the model and develop in the same manner while under the same environmental conditions. Growth equations of hydrometeor particles with similar radii are solved together and the particle size distribution (PSD) f (r), which is the number density sorted by the radius of the hydrometeors, is then predicted. Numerous observational studies have shown that liquid hydrometeor particles have three major modes for the PSD (Pruppacher and Klett [1997]): particles in the condensational growth mode (which have a mode radius around 10 μm), particles in the collisional growth mode (which have a mode radius larger than 100 μm), and the drizzling mode (which is a transitional mode between the two aforementioned modes).

There are two main methods that can be used to predict the time evolution of the PSD. One is the multi-moment bulk method, which predicts the integrand of each mode by approximating the modes via the gamma distribution (e.g., Milbrandt and Yau [2005b]; Seifert and Beheng [2006]):

$$f\left(r\right)=K{r}^{v}exp\left(-\lambda r^{\mu }\right)\text{,}$$

where κ, λ, μ, and ν are the diagnosed parameters. The k-th moment of the PSD, M^(k), is defined as

$$M^{\left(k\right)}\equiv {\displaystyle {\int }_0^{\infty }f\left(r\right)r^k\mathit{dr}\text{,}}$$

and the time evolution of the prognostic moments is evaluated by integrating the growth rates of a single hydrometeor particle over the entire range of the PSD (e.g., Seifert and Beheng [2006]). The other is the spectral bin method, which solves growth equations by discretizing the PSD with tens of bins covering these modes (Khain et al. [2000], [2008]). Because of their computational efficiency, the multi-moment bulk methods are widely used for 3D non-hydrostatic simulations.

The multi-moment bulk method is categorized by the number of its prognostic moments: the single-moment bulk method predicts the mass concentration of hydrometeors (Kessler [1969]; Lin et al. [1983]; Rutledge and Hobbs [1983]; Walko et al. [1995]; Grabowski [1998]; Hong et al. [2004]; Thompson et al. [2008]; Tomita [2008b]), and the double-moment bulk method predicts the number concentration in addition to the mass concentration (Meyers et al. [1997]; Feingold et al. [1998]; Morrison et al. [2005]; Seifert and Beheng [2006]; Phillips et al. [2007]; Morrison and Gettelman [2008]; Lim and Hong [2010]).

The four parameters used to represent the PSD (κ, λ, μ, and ν) are diagnosed by prognostic moments. Hence, as the number of the prognostic moments is increased, a variety of PSD shapes can be represented. For example, increases in the mode radius and skewness of the PSD are observed in collisional growth (Berry and Reinhardt [1974]). These distortions of the PSD are important for initiation of precipitation (Seifert and Beheng [2001]). Furthermore, once precipitation occurs, the gravitational sedimentation of the hydrometeor particles induces further distortion of the PSD due to the gravitational size sorting mechanism and the time evolution of the PSD affects the subsequent chain of particle growth (Wacker and Seifert [2001]; Milbrandt and Yau [2005a]; Milbrandt and McTaggart-Cowan [2010]). In addition, higher order moments are required for the closure of the governing equations in the collisional processes (e.g., Seifert and Beheng [2001]). Thus, an increase in the number of prognostic moments is linked to improvements in the time evolution of the PSD and complexity of the cloud microphysics modeling. Additionally, complexity and calculation costs increase as the number of prognostic moments increases.

For climate studies, collisional processes are particularly important to the radiative budget because such processes dominate the persistence of warm-phase clouds due to the rapid growth rate of cloud droplets during the precipitating process, along with the strong cloud albedo perturbations caused by the spatial and temporal variabilities of their growth rates (Albrecht [1989]; Lohmann and Feichter [2005]). The double-moment bulk method enables us to explicitly evaluate the dependence of the collisional growth rate on the mode radii (Berry and Reinhardt [1974]; Feingold et al. [1998]; Khairoutdinov and Kogan [2000]; Seifert and Beheng [2001]). In addition, cloud radiative forcing can be more accurately estimated because the cloud optical properties are calculated using the mass concentration and number concentration (Hansen and Travis [1974]). In contrast, the single-moment bulk method assumes that the radii of the hydrometeors are empirical values. Empirical values do not reproduce the temporal and spatial variability of atmospheric states, and their validity has not been confirmed for various types of climate change, such as global warming conditions. The triple-moment bulk method, which was recently developed and used (Milbrandt and Yau [2005a], [b]; Seifert [2008]; Shipway and Hill [2012]), predicts or diagnoses the radar reflectivity factor in addition to the aforementioned two moments and is applied to the regional simulations of severe storms (Milbrandt and Yau [2005a], [2006]). Of course, the higher order moment bulk methods also contain empirical formulations and parameters that can be used to calculate cloud growth and optical properties. However, previous process studies have shown that the differences of simulated results between double-moment and triple-moment bulk methods are smaller than those between single-moment and double-moment bulk methods (e.g., Milbrandt and Yau [2005a]; Shipway and Hill [2012]). In addition, it has also been demonstrated that a well-optimized double-moment method can capture the characteristics of cloud growth simulated by a spectral bin method (Seifert and Beheng [2001]; Seifert et al. 2006; Seifert [2008]; Milbrandt and Yau [2005a], [2006]). These facts indicate that the double-moment bulk method is practical and applicable in terms of cost performance of numerical simulations.

Cloud microphysics schemes in NICAM

NICAM initially adopted a single-moment bulk cloud microphysics scheme that was proposed by Grabowski ([1998]) in order to explore the multi-scale interaction of moist convective systems because the organizations of deep convective systems do not depend on the details of cloud microphysics schemes. Instead, meso-scale circulations are reproduced if the latent heat release in updrafts and the evaporation cooling caused by rain drops are considered by predicting the mass concentration of hydrometeor particles, as introduced in the single-moment bulk scheme. In the latest version, the schemes proposed by Kessler ([1969]), Lin et al. ([1983]), Hong et al. ([2004]), and Tomita ([2008b]) are implemented. In particular, NICAM single-moment scheme with six water categories proposed by Tomita ([2008b]) (hereinafter, referred to as NSW6) has been well evaluated by comparisons with satellite observations (Satoh et al. [2010]; Kodama et al. [2012]; Hashino et al. [2013]; Roh and Satoh [2014]) and, hence, is used by default. Because of the relatively low calculation cost, the use of single-moment schemes enables us to conduct very high-resolution global non-hydrostatic simulations at g-level 13 ($\overline{\Delta }$ = 870 m; Miyamoto et al. [2013]), or for longer integrations up to several decades.

Recently, Seiki and Nakajima ([2014]) implemented a double-moment bulk cloud microphysics scheme (NICAM double-moment scheme with six water categories; hereinafter, referred to as NDW6) into NICAM. Seiki et al. ([2014]) demonstrated that the NDW6 scheme well reproduces cloud radiative forcing by calculating cloud ice number concentration with using accurate databases of the cloud optical properties. Figure 6 shows the simulated cloud optical thickness (COT) and outgoing shortwave radiative flux at the top of the atmosphere that were obtained using the NDW6 scheme, along with the observed outgoing shortwave flux at the top of the atmosphere that was obtained from the satellite product of CERES SYN1deg_Ed3A (Wielicki et al. [1996]; Kato et al. [2011]). The simulation showed that convective clouds with a high COT were scattered over the tropics and formed regional cloud clusters (e.g., over the central Pacific). For example, Typhoon Fengshen remained over the Philippines and high COT clouds were distributed along its vortical structure. Furthermore, synoptic scale disturbances were organized over the mid-latitude regions and optically thick regions were located along the front. A shortwave flux corresponding to these high COT regions was largely reflected by clouds. These features agree well with the satellite observation. As shown by the NDW6 scheme, the impact of the cloud radiative forcing, which results in improvements to the climatological biases in NICAM, should be evaluated in future studies.

Figure 6

Simulated results for the 14 km mesh NICAM with the double-moment cloud microphysics scheme. Simulated results for the double-moment cloud microphysics scheme for the 14 km mesh NICAM on 20 June 2004: (a) simulated cloud optical thickness at a wavelength of 550 nm at 00:00 UTC, (b) daily averaged outgoing shortwave flux (W m⁻²) at the top of the atmosphere, and (c) the daily averaged outgoing shortwave flux (W m⁻²) at the top of atmosphere, as produced by the CERES SYN1deg_Ed3A product.

1.3.2 Subgrid-scale turbulence of the planetary boundary layer

Boundary layer turbulence plays an important role in driving large-scale circulations of the atmosphere as it transports heat upward from the Earth's surface in the form of sensible and latent heats, which are converted from incoming solar energy at the surface (Trenberth et al. [2009]). To simulate atmospheric circulations, a parameterization scheme that computes the effects of such turbulent transport in the subgrid scale of a numerical model needs to be incorporated, unless the resolution is high enough to resolve turbulence.

NICAM adopts the MYNN scheme, which is a modified version of the original Mellor-Yamada (MY) model (Mellor and Yamada [1982]; Nakanishi and Niino [2006], [2009]), for such turbulent transport models (Noda et al. [2010]). There are three major advantages to employing the MY-type model: 1) the MY-type model can consistently and simultaneously compute subgrid-scale condensation and resultant turbulent fluxes by assuming the probability density functions regarding the inhomogeneity of temperature and water on a subgrid scale, 2) the model can consistently evaluate the impacts of the differences of turbulent moments, and 3) the model can run with low computational cost because it omits minor terms in the turbulent transport equation system. Based on the degree of approximation from the original equations, level 1 to level 4 can be included in the MY-type model. The level-2, level-2.5, and level-3 MYNN version models are installed in NICAM, and the level-2 model is used as a standard PBL scheme.

Figure 7 compares the impacts of these different level models on the spatial distributions of simulated clouds. As can be seen in the figure, all turbulent models properly simulate the spatial characteristics of, for example, high and mid-level clouds, even though the magnitude of the high clouds tends to be larger than those seen in satellite observations. For low clouds, the level-2 model properly represents the characteristics of the regional contrast between a clear and cloudy sky, such as those seen over the eastern and western regions of the Pacific Ocean. Conversely, the low clouds that systematically develop along and off the Peruvian coast are not simulated well. Compared with the results of the level-2, level-2.5 and level-3 models produce a much smaller fraction of low clouds. It is speculated that these two models calculate more vertical turbulent heat transports of heat from the boundary layer to the free atmosphere, which results in a less humid boundary-layer state.

Figure 7

Comparison of the cloud amounts simulated by using different turbulent models to the satellite observation. The comparisons of the cloud amounts from the ISCCP satellite data, level-2, level-2.5, and level-3 MYNN turbulent schemes were obtained by using the 14 km mesh NICAM. The NICAM data shown are averages during 29 July 1982 and 15 August 1982, and that of the satellite observation are the climatology for July and August.

Due to the remarkable progress of computational technology, the model resolution of NICAM is steadily being increased, with the highest resolution to date having been achieved when a horizontal mesh size of 870 m was used (Miyamoto et al. [2013]). For a global non-hydrostatic model, boundary-layer scale turbulence becomes explicitly and partially calculatable around this resolution, but will not be sufficiently resolved until the mesh size is less than 100 m (Khairoutdinov et al. [2009]; Moeng et al. [2009]). It is a challenging issue to parameterize turbulent transport processes in such an intermediate resolution (Wyngaard [2004]), and continuous efforts are needed to improve the subgrid-scale turbulence parameterization schemes along with shallow cumulus convection.

1.3.3 Radiation

Atmospheric radiative transfer models calculate heating rates in the atmosphere and radiative flux at the Earth's surface using the 3D distribution of clouds/gases/aerosols and land/ocean surface properties with the atmospheric state. NICAM adopts the broadband radiative transfer model ‘MSTRN’, which is based on the discrete ordinate method with a delta two-stream approximation and the correlated k-distribution method. The MSTRN model calculates gas absorption (H₂O, CO₂, O₃, N₂O, CH₄, and O₂) and the scattering/absorption of hydrometeors and aerosols (dusts, sea salts, carbonaceous, and sulfates) over an arbitrary time interval and then outputs short- and longwave radiation fluxes. A basic description of this model is written in Nakajima et al. ([2000]). In the latest package, ‘MSTRN-X’ (Sekiguchi and Nakaima [2008]), the gas absorption processes are improved and successfully used to reduce the error of the heating rate. The scheme was introduced into NICAM through the AGCM version implemented for MIROC from the original one-dimensional (1D) model. The calculation of solar insolation (Berger [1978]) was added to the AGCM version. The other cloud overlapping method (maximum-random overlapping) was incorporated since this version of the radiation scheme was included in MIROC3.2 (Hasumi and Emori [2004]).

The calculation cost of the radiation scheme is very high in the AGCM, because both computing speed and precision are required for the scheme. The latest MSTRN-X scheme improved the optimization method in order to reduce the integration points in the correlated k-distribution approximation without decreasing accuracy. In the NICAM default setting, the spectrum between 0.2 and 200 μm is divided into 29 spectral bands with 111 integration points. Cloud-radiation feedback is one of the most important issues in the study of climate change. The high-resolution experiments performed by NICAM have an advantage in that they reduce the uncertainties of the horizontal/vertical cloud distributions. In the recent studies described above, both the cloud scheme and radiation scheme are improved by considering the detailed size distribution of the hydrometeors and non-sphericity of the ice particles. This cooperative development of the cloud and radiation schemes will facilitate improvements in the future development of aerosol and chemistry modules (see ‘Aerosol and chemistry modules’ section).

NICAM also introduces the ISCCP simulator (Klein and Jakob [1999]; Webb et al. [2001]), with which the distribution of every cloud category can be calculated online. The results from the simulator have been used in numerous studies (Iga et al. [2007a]; Suzuki et al. [2008]; Collins and Satoh [2009]; Noda et al. [2010]; Sohn et al. [2010]; Satoh et al. [2012b]; Miyamoto et al. [2013]; Tsushima et al. [2014]). In some studies, new satellite simulator packages have been applied offline to the NICAM dataset. Kodama et al. ([2012]) used a series of satellite simulators that were part of the Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulator Package (COSP) (Bodas-Salcedo et al. [2011]), and Hashino et al. ([2013]) used the Joint Simulator for Satellite Sensors (J-simulator). Therefore, the online use of these packages should also be taken into consideration.

1.3.4 Aerosol and chemistry modules

Currently, NICAM is implemented with a global 3D aerosol transport-radiation model, the Spectral Radiation-Transport Model for Aerosol Species (SPRINTARS) (described fully in Takemura et al. [2000], [2002], [2005], [2009]; Goto et al. [2011a]). The SPRINTARS module predicts the mass mixing ratios of the main tropospheric aerosol species, i.e., carbonaceous (black carbon and organic matter), sulfate, soil dust, sea salt, and the precursor gases of sulfate (sulfur dioxide and dimethylsulfide (DMS)). The aerosol transport processes include emission, advection, diffusion, sulfur chemistry, wet deposition, dry deposition, and gravitational settling. Emissions of soil dust, sea salt, and DMS are calculated using the internal state of the model, whereas the emissions for the other primary aerosols and the distributions of oxidants (hydroxyl radical, ozone, and hydrogen peroxide) for the sulfur chemistry are prescribed by externally assumed data. In the module itself, meteorological fields such as clouds and precipitation, which are critical to determining aerosol lifetimes, are provided from the outer module after the cloud microphysics calculations are computed at each time step. In contrast, the calculated mass mixing ratios of the aerosols in each grid cell are calculated by the aerosol module SPRINTARS and then returned back to the outer module. In addition, aerosol variables are forwarded to the radiative transfer module, MSTRN-X (Sekiguchi and Nakaima [2008]), in order to represent the direct aerosol effects, and also to cloud schemes such as NSW6, to represent the indirect aerosol effects.

SPRINTARS, which was originally coupled to MIROC (Watanabe et al. [2010]) with a low spatial resolution of 100 to 300 km, has compared well with various measurements (Goto et al. [2011a], [b], [2012]) and other global aerosol models participating in the AeroCom international project (e.g., Kinne et al. [2006]). The results of the SPRINTARS module coupled to NICAM were also shown to generally agree with the satellite observations of aerosols and clouds for versions based on low spatial resolution of g-level 5 ($\overline{\Delta }$ = 220 km) (Dai et al. [2014a], [b]; Goto [2014]) and high spatial resolution of g-level 10 ($\overline{\Delta }$ = 7 km) (Suzuki et al. [2008]). Using Earth Simulator, Suzuki et al. ([2008]) conducted a global simulation with the finest resolution to date and showed the global aerosol and cloud distributions for July 2006 (Figure 8). Most of the aerosol plumes were in generally good agreement with the satellite observations, except for North America, Australia, and most remote oceans. Especially, over North America, the values of NICAM-simulated aerosol optical thickness were underestimated, possibly because of insufficient emissions from biomass burning and anthropogenic sources, or the lack of secondary organic aerosols formed from isoprene in forests (Henze and Seinfeld [2006]). Through the aerosol-cloud interaction, the cloud droplet effective radius generally decreases, as the aerosol optical thickness increases, which were generally captured by both NICAM and MODIS. The cloud droplet effective radii simulated by NICAM were smaller than those obtained from MODIS, probably because those retrieved from 2.2 μm (shown in Figure 8) correspond to those captured below the cloud top heights where the cloud droplet effective radii are larger than those at cloud top height (Nakajima et al. [2010]). Very recently, regional aerosol simulations that used the stretched grid system with a high spatial resolution of approximately 10 km were performed for the Kanto region of Japan under the MEXT/RECCA/SALSA project (SALSA [2014]; Goto et al. [2014]). These simulations were carried out at relatively low computational cost by using two supercomputers at the University of Tokyo (SR16000 and FX10).

Figure 8

Distributions of aerosol optical thickness and cloud droplet effective radius from the NICAM-SPRINTARS simulations. Global geographical distributions of (a, c) aerosol optical thickness and (b, d) cloud droplet effective radius from (c, d) the NICAM-SPRINTARS simulations in comparison to those obtained from (a, b) the MODIS satellite observations for 1 to 8 July 2006 (cited from Suzuki et al. [2008]). The unit of cloud droplet effective radius is micrometers.

The CHASER atmospheric chemistry module has also been implemented into NICAM (SALSA [2014]). This module primary focuses on tropospheric chemistry while considering the chemical cycle of O _x -NO _x -HO _x -CH₄-CO relation to the oxidation of volatile organic compounds (VOCs) and halogen chemistry in order to calculate the concentrations of 92 chemical species with 262 chemical reactions (58 photolytic, 183 kinetic, and 21 heterogeneous reactions) (Sudo et al. [2002a], [b]; Watanabe et al. [2011]). Using the CHASER module, the mixing ratios of atmospheric gases, which mainly include ozone, NO _x , CO, SO₂, VOCs, and sulfate aerosols, can be predicted. The transport processes of the chemical tracers include emission, advection, diffusion, chemical reaction, wet deposition, and dry deposition. The emissions of lighting NO _x and DMS are calculated by using the internal parameters of the model, whereas the emissions of the other gases are prescribed. The photolysis rates are calculated online by using the temperature and radiation fluxes computed in the radiative transfer module. Like the aerosol module, the meteorological field information is input from the outer module after the cloud microphysics calculations and to the chemistry module at each time step. In contrast, the calculated mass mixing ratios of the tracers in each grid cell within the chemistry module are returned back to the outer module. In addition, ozone, N₂O, and CH₄ are calculated online and used in the radiative transfer module to represent the atmospheric warming effects.

The distributions of the trace gases obtained by CHASER, which is currently coupled to MIROC at a low spatial resolution of approximately 300 km, are generally compared to the observations (Sudo et al. [2002b]; Sudo and Akimoto [2007]; Nagashima et al. [2010]) and to other chemistry models under international model comparison projects such as the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP) (e.g., Lamarque et al. [2013]). Under the MEXT/RECCA/SALSA project, the results of CHASER, along with the globally uniform and regionally stretched grid systems, were compared with the observations (SALSA [2014]). Furthermore, the chemical module has been currently expanded to a fully coupled aerosol-chemistry module by coupling CHASER with SPRINTARS, as developed in Watanabe et al. ([2011]) and used by MIROC.

The unified module, consisting of the fully coupled aerosol-chemistry, coupled to NICAM (hereinafter referred to as NICAM-Chem) can explicitly consider various interactions between aerosols and gases, as shown in Figure 9. For example, as an aspect of the aerosol simulation, the unified module can precisely calculate the sulfate formation with online oxidants such as hydroxyl radicals and ozone, whereas SPRINTARS uses assumed oxidants. In contrast, as an aspect of short-lived gases, calculations of the heterogeneous reactions onto the aerosol surface are expected to provide the greatest advantage when using the unified model. In addition to the sophisticated interactions between aerosols and gases, the chemistry of nitrate will be implemented into NICAM-Chem by using a thermodynamic equilibrium condition and secondary organic aerosols (SOAs). As a result, the new unified chemistry module is expected to provide better distributions of aerosols and gases, and, in the near future, we expect to be able to provide such distributions globally with a high spatial resolution of less than 10 km.

Figure 9

Relationship between the aerosol, gas, cloud, precipitation, and radiation fields for NICAM-Chem. Relationship between the aerosol, gas, cloud, precipitation, and radiation fields for the unified aerosol-chemistry module (NICAM-Chem). The aerosol and gas modules correspond to SPRINTARS and CHASER, respectively.

1.3.5 Land model

We implemented two land surface schemes in NICAM in order to calculate the water and energy balance over land. The first scheme is a simple bucket land surface model (Kondo [1993]), and the other involves the minimal advanced treatments of surface interaction and runoff (MATSIRO) model (Takata et al. [2003]). When used with NICAM, the simple bucket land surface model, which was first used by Miura et al. ([2007a]), has three vertical layers for predicting soil temperature and one layer for predicting the soil moisture.

The MATSIRO scheme is the same as that originally employed in MIROC (Hasumi and Emori [2004]). Kodama et al. ([2012]) were the first to use MATSIRO with NICAM. This scheme has five vertical layers for predicting soil temperature and soil moisture, as well as a snow scheme with thee vertical layers.

Even in a several-kilometer grid simulation, land use and land cover cannot be regarded as homogeneous in one grid cell if it is to be used to evaluate the land surface heterogeneities smaller than the horizontal grid scale. Accordingly, we are currently implementing a mosaic treatment of the land surface models in order to improve the accuracy of the water and energy fluxes at the land surface. We will also attempt to incorporate an urban canopy model in order to simulate the global urban energy balance.

1.4 Assimilation

1.4.1 Local ensemble transform Kalman filter applied to NICAM: NICAM-LETKF

The introduction of a data assimilation system into NICAM is a challenge because it might enable us to produce a high-resolution dataset with O (km) mesh size that can be used for improved better numerical simulations and research analysis. Data assimilation is an analysis technique that provides an accurate estimate of the atmospheric state by using observation data and model forecasts. It plays an important role in improving numerical weather prediction by providing an accurate initial condition. The ensemble Kalman filter (EnKF) is an advanced data assimilation method that uses flow-dependent forecast error statistics that are represented by the ensemble prediction (Evensen [1994]). For a linear system, the Kalman filter (KF) gives a minimum mean square error if the error is a Gaussian distribution (Kalman [1960]). However, it is difficult to apply the KF because the computational cost for an atmospheric model with a high degree of freedom (order of 10⁷) is high. To resolve this problem, Evensen ([1994]) proposed the EnKF, which approximates the KF by using the forecast spread of the ensemble prediction. Furthermore, the local ensemble transform Kalman filter (LETKF) (Hunt et al. [2007]) was developed based on the ensemble transform Kalman filter (ETKF) (Bishop et al. [2001]) and used an algorithm that was designed for suitability with parallel computers by taking advantage of the independent local analyses of the local ensemble Kalman filter (LEKF) (Ott et al. [2004]). For details on the mathematical formulation and implementation of these approaches, refer to each reference.

Previous studies have applied the LETKF to regional and global atmospheric models and have shown promising results (Miyoshi and Aranami [2006]; Miyoshi and Yamane [2007]). Kondo and Tanaka ([2009]) applied the LETKF to NICAM in order to develop NICAM-LETKF, which was the first LETKF to be applied to a global non-hydrostatic model. They then investigated the feasibility and stability of NICAM-LETKF using a perfect model and determined that the system works appropriately for a realistic non-hydrostatic model (NICAM). However, the horizontal resolution used in their study was not very high (g-level 5 or $\overline{\Delta }$ = 220 km), where the model would behave hydrostatically. Additionally, it is unclear whether NICAM-LETKF will work stably when real observations are used. Therefore, a future task will be to investigate data assimilation with a high-resolution global non-hydrostatic regime using real observations. This would not only contribute to an improvement in weather prediction but also to the research on multi-scale interaction, by using the product of the assimilation. We are now developing a version of NICAM-LETKF that will be based on the LETKF code that was developed and continuously improved by Miyoshi (e.g., Miyoshi and Kunii [2011]) and are currently testing assimilation of precipitation data observed by satellites in order to obtain a fine-mesh precipitation database using the method proposed by Lien et al. ([2013]).

1.4.2 NICAM-based transport model: NICAM-TM

An accurate estimation of the past and current CO₂ budget for the Earth's surface is needed for reliable predictions of carbon cycle changes under global warming conditions. Therefore, an inversion method coupled with a tracer transport model was used to estimate the spatial and temporal variations of CO₂ sources and sinks based on observations of CO₂ concentrations (e.g., Tans et al. [1990]). In an inversion analysis, we quantitatively estimated regional CO₂ fluxes by considering the a priori estimate and uncertainties of the CO₂ fluxes, which were derived from our present understanding of CO₂ source/sink mechanisms and uncertainties of observations (including model representative errors). This analysis is based on Bayesian theory and is similar to data assimilations for meteorological forecasts (Tarantola [2005]). However, in contrast to the observations of meteorological parameters, CO₂ concentrations are not sufficiently observed around the globe. Therefore, estimates of CO₂ sources and sinks are not well constrained by actual observations and are extremely sensitive to model transport errors (e.g., Gurney et al. [2002]).

In response to the need for an accurate transport model, Niwa ([2010]) developed the NICAM-based transport model (NICAM-TM). For transport simulations of CO₂, NICAM has an advantageous property: consistency with continuity (CWC) (Jöckel et al. [2001]; Gross et al. [2002]). Under the CWC condition, the mass conservation and Lagrangian conservation of a tracer are simultaneously guaranteed; that is, the volume integral of mass is conserved and a constant specific mass of a tracer is maintained along flow trajectories. The ability to achieve the CWC property is attributed to the fact that NICAM uses the finite volume method for its meteorological field integration and performs tracer transport on a common dynamical frame (Niwa et al. [2011a]). The tracer mass conservation is strictly required for analyses of the CO₂ budget, which is the main target of inversion studies. Furthermore, the Lagrangian conservation property is imperative for CO₂ because it is an abundantly existing tracer (the global average CO₂ concentration was about 390 ppm for 2013) and the concentration changes that are analyzed are very small (no more than a few tens of ppm).

NICAM-TM consists of online and offline transport models, as well as an adjoint model for the tracer transport. In the online model, tracer transport is calculated concomitantly but independently with the integration of the meteorological field. More specifically, the tracers other than those for water are passive. Compared to the online model, the offline model is computationally inexpensive because it only calculates the tracer transport by using meteorological data calculated by the online NICAM-TM and stored in a disk beforehand. The adjoint model was developed based on the offline model, which calculates the sensitivities to certain tracer variables backward in time. Because long-term and multiple simulations are required for the inverse calculation, NICAM-TM is currently run with a relatively low-resolution for CO₂ transport simulations; for the most part, a resolution of g-level 5 or 6 is used. However, after conducting intercomparison experiments (Law et al. [2008]; Patra et al. [2008]; Niwa et al. [2011b]), we found that the performances of CO₂ transport obtained are comparable to those of other models. Generally, the modeled winds are nudged towards the analyzed data (e.g., JMA Climate Data Assimilation System (JCDAS); Onogi et al. [2007]) in order to match the simulated CO₂ concentrations to the observed concentrations.

Niwa et al. ([2012]) performed an inversion analysis by using the synthesis inversion method (Enting [2002]) with NICAM-TM. This was the first inversion study in which aircraft data were extensively used to constrain the global and regional CO₂ budgets. The aircraft data were obtained from the Comprehensive Observation Network for Trace Gases by Airliner project (CONTRAIL) (Machida et al. [2008]; Matsueda et al. [2008]; Sawa et al. [2012]). By adding the CONTRAIL data (Figure 10a) to the conventionally used surface observation network, Niwa et al. ([2012]) successfully reduced uncertainties of the CO₂ flux estimates for Asia. Moreover, this study revealed the existence of a strong summer uptake by the biosphere in South Asia (Figure 10b), where the CO₂ observations from CONTRAIL in the middle to upper troposphere have a predominant constraint on the flux estimation. Because biosphere models have not shown such a strong summer uptake (Patra et al. [2013]), Niwa et al. ([2012]) suggested the existence of some unknown mechanisms of the CO₂ sources and sinks in this region. To further improve the CO₂ inversion, a four-dimensional (4D) variational method system is now being developed by using the adjoint model of NICAM-TM. Using this system, a large amount of CO₂ concentration data that have been recently available from regular aircraft observations and satellite observations (e.g., Greenhouse Gases Observing Satellite (GOSAT)) will be exploited to further constrain the estimates of CO₂ sources and sinks.

Figure 10

Simulated CO ₂ concentration field and monthly variation of the CO ₂ flux by NICAM-TM. (a) Mean CO₂ concentration field at 250 hPa for June 2007, as simulated by NICAM-TM, and the locations of the CONTRAIL measurements for 2007 (white dots). (b) Monthly variation of the CO₂ flux averaged for 2006 to 2008, estimated by the inversion analysis of Niwa et al. ([2012]) (this figure is slightly modified from Figure 5d in Niwa et al. [2012]). The blue and red lines are the CO₂ fluxes (positive values indicate sources and negative values indicate sinks) estimated by only using surface data and by using both surface and CONTRAIL data, respectively, and the colored shades and error bars denote the range of the estimated flux uncertainty and standard deviation for 2006 to 2008, respectively.

1.5 Models related to NICAM

1.5.1 Stretch NICAM and diamond NICAM

As described in ‘Grid configuration and advection scheme’ section, NICAM was originally designed for use with an icosahedral grid, which covers a sphere of the Earth quasi-uniformly and has a grid system based on 10 large diamonds. As extensions to the original configuration, NICAM can also be applied to a non-uniform resolution distribution grid or a regional grid by changing the positions of the grid points, the total number of diamonds, or the topology of their combination.

NICAM can also be adapted to simulations that target a particular area. Two approaches to performing a partially high-resolution simulation are typically used with NICAM. One includes the stretched icosahedral grid system developed by Tomita ([2008a]). In this approach, the icosahedral grid is stretched using the Schmidt transformation method and the grid size becomes gradually finer as grid points are concentrated onto the area of interest. In the stretch NICAM, the inhomogeneity of the horizontal grid sizes requires resolution-dependent parameters for physical schemes if the global domain is of interest. However, if a limited region is of interest, such as within regions with finer mesh sizes, resolution-dependent parameters are not required. The stretch NICAM is widely used for studies on tropical convective systems (Satoh et al. [2010]; Satoh and Kitao [2013]; Roh and Satoh [2014]), tropical cyclones (Yanase et al. [2010a], [b]; Satoh et al. [2013]; Arakane et al. [2014]), and MJOs (Nasuno [2013]).

The other approach is similar to the use of a regional model. Only 1 of the 10 diamond regions that make up the icosahedral grid is used for the simulation (hereinafter, the diamond NICAM). The simulation domain of the diamond NICAM consists of 1 diamond (2 triangles of an icosahedron), even though the original NICAM global domain consists of 10 diamonds (20 triangles). The diamond NICAM is used as a non-hydrostatic regional model. By using this approach, we can share the same code for both the global domain and limited-area models. This is advantageous for the maintenance of the source code because all components of NICAM are used in common in both models.

To enable the limited-area simulations using the diamond NICAM, we developed a nudging scheme in which the strength of the nudge depends on the distance from a certain point on the globe. Figure 11 shows the simulation domain of a regional climate simulation targeted on Japan. The lateral boundary of the simulation domain (outer edge of the gray circle in Figure 11) is forced by six-hourly reanalysis data. Verification of the diamond NICAM is still under investigation and the results of test simulations performed using the diamond NICAM were compared with the results of the high-resolution experiments using the stretch and global NICAM. The comparison results showed good agreement with those produced by the original global NICAM simulation. The temporal evolution of a convective system, such as the locations of triggering, development, and decay, was also well simulated during the integration.

Figure 11

Simulation domain of the diamond NICAM targeted on Japan. Light blue dots indicate the grid points. Reanalysis nudging was applied to the outer side of the gray circle. This grid is based on g-level 5. Stretching parameter (β in Equation 8, Tomita et al. [2008]) of grid transformation is 10. The minimum grid interval is about 59 km.

1.5.2 Plane NICAM

Although the stretch and diamond NICAM can be used for limited-area simulations, an f-plane system is useful when conducting some types of numerical simulations, such as idealized simulations on an f-plane and simulations with a narrow region in which the curvature of the Earth is not important. For this reason, an f-plane model (hereinafter, the plane NICAM) whose dynamical core is the same as the original NICAM was developed. To save maintenance costs, the plane NICAM used the same code as the global NICAM, which also makes knowledge gained by the plane NICAM immediately and directly applicable to the global NICAM.

In the plane NICAM, the shape of the calculation domain is a diamond composed of two regular triangles with a double periodic lateral boundary condition. Figure 12 shows the configuration of the calculation domain for the plane NICAM, as indicated by black bold lines. The outward fluxes across side AB are equivalent to the inward fluxes across side DC. The red, green, blue, and yellow areas inside the diamond area surrounded by the black bold lines correspond to the hatched areas of the same respective color. Therefore, if the double periodic boundary condition is applied to the diamond-shaped domain with side length L, the calculation domain is equivalent to the regular hexagonal domain with side length l = L/√3, as represented by the hatched area in Figure 12. Regular hexagonal geometry is better suited than square or regular triangle geometry in terms of isotropy.

Figure 12

A schematic image of the calculation domain for the plane NICAM. A schematic image showing the relationship between the calculation domain and a unit cell for the plane NICAM. The diamond-shaped area (ABCD) surrounded by black bold lines indicates the calculation domain, and the hatched regular hexagonal area indicates the unit cell.

As an example of the use of the plane NICAM, we briefly show a result of an idealized tropical cyclone simulation on an f-plane (Ohno and Satoh [2014]). The model was initialized with an axisymmetric cyclonic vortex that had a maximum azimuthal wind of 12 m s⁻¹ as formulated by Rotunno and Emanuel ([1987]). The initial thermodynamic structure of the unperturbed model atmosphere is defined by moist adiabatic lapse rate up to 17 km and a constant temperature above this height. The length of the sides of the numerical domain is L = 3,000 km, and a horizontal resolution of 2.7 km is used for the entire domain. We used 50 vertical levels up to 45 km where the intervals of the vertical levels become smaller in the boundary layer. The SST is fixed at 31°C and an f-plane with a constant Coriolis parameter at 18 N is assumed. After an initial shock, the minimum sea level pressure continued to drop for 140 h and a quasi-steady stage was achieved. The cyclone reached a minimum central surface pressure of 931 hPa at about 160 h. A plan view of the distribution of the hourly precipitation amounts after 128 h of the time integration is shown in Figure 13, which indicates that the tropical cyclone has an eye, an eyewall, inner core rainbands, and distant rainbands. The array of rainbands and the eyewall are typical features of a tropical cyclone (Houze [2010]).

Figure 13

A horizontal distribution of precipitation for a tropical cyclone simulation by the plane NICAM. A horizontal distribution of hourly precipitation amount by the plane NICAM after a time integration of 128 h for a 2.7 km mesh tropical cyclone simulation.

1.5.3 Further modified grids

In addition to the original icosahedral grid, alternative spherical grids denoted by extended triangular meshes (XTMS) (Iga [2014]a) are also implemented on NICAM. The mesh topologies of these grids are different from the topology of an icosahedral grid and are generated by three processes: grid relaxation via spring dynamics, transformation via an analytical function, and the Schmidt transformation. By changing mesh topologies and transformation functions, it is possible to obtain a grid with non-uniform distributions of resolution on a sphere, which is applicable to various situations. For example, the resolution can be increased at a particular region or along a tropical zone as shown by Figure 14a,b. The grid structure shown by Figure 14a looks similar to that of the stretch NICAM, but the variability of the horizontal mesh spacing is smaller and the variability of the resolution is smoother than that of the stretch NICAM. The resolution is approximately proportional to the combination of map factors of two polar stereographic projections. As for the grid structure shown by Figure 14b, resolution is enhanced along the equatorial belt and is approximately proportional to the combination of map factors of two Lambert conformal conic projections and one Mercator projection (hereinafter, LML grid). The LML grid is more useful for studying the multi-scale structure of convective systems in the tropics. The performance of the LML grid was examined by aqua planet experiments (Iga [2014]b).

Figure 14

Grid configurations of the regionally enhanced grid and the tropically enhanced grid. Grid configurations of (a) the regionally enhanced grid and (b) the tropically enhanced grid (after Iga [2014]a).

1.5.4 Coupler

Currently, high-resolution NICAM is mainly used for short-term simulations such as for tropical cyclones and intra-seasonal variability. Since deep sea circulations are not important for these simulations, a simple mixed-layer ocean model is implemented in NICAM. However, in order to cover wider areas of research, such as climate projection, coupling with an ocean model is necessary when NICAM is used as an atmospheric component of an Earth system model. To couple an atmosphere model and an ocean model, we used a coupler called Jcup (Figure 15). Because Jcup connects models with different grid systems, it can be used as a grid transformation tool. Thus, in addition to ocean coupling, Jcup is used as an I/O tool of NICAM by converting the output data from the icosahedral grid to the lat-lon grid which is generally used for analysis. In this section, the Jcup coupling library that was used as a core of the coupled model system will be introduced first, after which the ocean and I/O coupling will be described.

Figure 15

A schematic of the coupling system.

The Jcup coupling library is a collaborative work of JMA/MRI, JAMSTEC, and the Research Organization for Information Science and Technology (RIST) (Arakawa and Yoshimura [2009]; Arakawa et al. [2011]), that is freely available as an open source program at https://sites.google.com/a/rist.jp/jcup/. Jcup inherited the development experience and design concept of Scup, which was developed by Yoshimura and Yukimoto ([2008]). The most remarkable feature of Jcup is its wide applicability. Generally speaking, the applicability range of a coupler is restricted by the specific grid structures and interpolation schemes it supports. For example, OASIS version 4 supported logically rectangular (i.e., 2D structured) or reduced Gaussian grids, along with bilinear, trilinear, bicubic, nearest-neighbor, and 2D conservative interpolations (Redler et al. [2009]). This indicates that the models that utilize other grid systems, such as non-structured or substructured grids, cannot be coupled.

The spectral method has been widely used for global atmospheric models. However, this method is not suitable for massively parallel processor (MPP) supercomputers because all-to-all communication is required. This is one reason why substructured grid models such as NICAM were developed energetically by the research organizations of numerous countries. Therefore, it is also necessary to develop a coupler that can handle such models. Jcup is specifically designed to couple different grid systems, so users can (and must) implement their own interpolation code and prepare a ‘mapping table’ in advance. See Arakawa et al. ([2011]) for details.

Next, we will describe NICAM and ocean model coupling. Currently, the Center for Climate System Research (CCSR) Ocean Component Model (COCO) (Hasumi [2000], [2006]) is used for the ocean model coupled with NICAM. In the coupler of the NICAM-COCO coupled model, Jcup is used as a core library and is implemented via a package called Joint-shell, which contains definition of the grid structures, exchange data, and the interpolation code (Figure 15). The interpolation code is based on the first-order conservative remapping scheme in Jones ([1999]). The ‘mapping table’ that must be set up as external files for beforehand defines grid-to-grid relations and Jcup's weight for grid systems of coupled models. Furthermore, we also implement a utility program to calculate such table files as a part of the Joint-shell. Note that COCO adopts a tri-polar grid, in which the northern polar region grid points do not follow latitude-longitude lines. As a result, coupling NICAM and COCO requires coupling two substructured grid models, which is why Jcup is used as a core of the coupler.

To couple NICAM, a new module has been implemented to replace the mixed-layer ocean module that has already implemented. Most of the other parts, such as the dynamical core, cloud microphysics, and land models, are untouched. The physical quantities exchanged between NICAM and COCO are listed in Table 3. The most significant issue was the inconsistency of ocean definition points between the two models. More specifically, some points that are recognized as ocean points by NICAM are recognized as land points by COCO. At such grid points, no data are provided to NICAM. Inconsistencies such as this occur if the land-ocean distribution has been set independently between the two models. In the future, we plan to improve pre-process utilities in order to create a NICAM land-ocean map that will maintain consistency with the COCO map. However, to avoid this problem in the interim, we applied the mixed-layer ocean module, which is pre-implemented in NICAM to such grid points.

Table 3 Physical quantities exchanged between NICAM and COCO through the coupler

The icosahedral grid employed by NICAM cannot be used to analyze the results. For example, the calculation of the zonal mean values is not straightforward and the visualization tools assume a lat-lon grid in many cases. Because of this, we generally use a conversion tool called ‘ico2ll’ to move data from the icosahedral grid to the lat-lon grid for post-processing of simulations. However, in many cases, running ico2ll as a post-process results in a bottle neck because it requires numerous I/O transfers and cannot parallelize itself since it is necessary to create a single file for the lat-lon grid of each variable. Jcup allows NICAM-I/O coupling, which leads to an efficient I/O conversion, to be successfully achieved. More specifically, the I/O module converts the icosahedral grid to the lat-lon grid and is executed in parallel with NICAM. The implemented conversion schemes include a bi-linear interpolation, a control volume weighted average, and the nearest-neighbor method. A bi-linear interpolation is originally implemented method in the stand-alone version of ico2ll’, and this method has first-order accuracy. The control volume weighted average scheme is the same with the one used in NICAM-COCO coupling and also has first-order accuracy, but the global averaged fluxes are conserved. The nearest-neighbor method has less accuracy and conservation properties than the other two methods. This method was implemented due to the demands from NICAM users and is used to provide a quick look at the raw values by using the general latitude-longitude graphic tools, such as the Grid Analysis and Display System (GrADS), without an averaging process. These schemes can be set for each output data type and can be switched by altering the configuration file.

Figure 15 is a schematic of the coupling system described above. The coupling system is designed so that NICAM automatically detects the coupling pattern at runtime without taking any other configuration. For example, when COCO is executed in parallel with NICAM, subroutines for NICAM-COCO coupling are used, and if not, subroutines of the mixed-layer ocean model are called. I/O is also the same as the case of COCO, in which NICAM automatically sends output data to I/O only when the I/O component is executed.

2 Conclusions

This paper reviewed recent activities of the global non-hydrostatic model, NICAM, and described its development and the details of each component of the dynamics and physics, the strategy of assimilation, and the related models. In the ‘Introduction’ section, a review of other global non-hydrostatic models, scientific overviews of NICAM, and the current design/structure and a future timeframe of NICAM were described. In the ‘Dynamics’ section, horizontal grid issues and vertical resolution issues were discussed. In the ‘Physics’ section, components of cloud microphysics, subgrid-scale turbulence, radiation, aerosol and chemistry, and land surface were described. As for assimilation, NICAM-LETKF was presented, and NICAM-TM was introduced to describe the assimilation of the transport model. The ‘Models related to NICAM’ section showed a variety of ways that NICAM can be used, specifically, stretch NICAM, diamond NICAM, plane NICAM, the tropically enhanced grid model, and coupling with the ocean model or the I/O module.

A number of comprehensive NICAM-related review articles have been published (Satoh et al. [2008]; Satoh [2013]) since the 3.5 km mesh global non-hydrostatic experiments were conducted by Tomita et al. ([2005]) and Miura et al. ([2007b]). However, these were primarily limited to descriptions of the governing equations, and because many of the model components continue to evolve and develop as the entire package expands, we were requested to provide a comprehensive overview of NICAM activities and future projections. NICAM has been a pioneering global non-hydrostatic model for almost a decade, and now that other global non-hydrostatic models are emerging, we believe that our experiences and considerations will be useful for other high-resolution modeling groups. Mutually beneficial information exchanges between different modeling groups will be required for future fruitful studies during the forthcoming era of global non-hydrostatic modeling.