The non-hydrostatic global atmospheric model for CMIP6 HighResMIP simulations (NICAM16-S): Experimental design, model description, and sensitivity experiments

NICAM, a nonhydrostatic global atmospheric model with an icosahedral grid system, has been developed for nearly two decades. This paper describes NICAM16-S, the latest stable version of NICAM (NICAM.16) modified for 15 Coupled Model Intercomparison Project Phase 6 (CMIP6). Major updates from NICAM.12, a previous version used for climate simulations, include updates of a cloud microphysics scheme and a land model, an introduction of natural and anthropogenic aerosols and a subgrid-scale orographic gravity wave drag, and improvement of coupling between cloud microphysics and radiation schemes. External forcings were updated to follow a protocol of CMIP6 High Resolution Model Intercomparison Project (HighResMIP). A series of short-term sensitivity experiments were performed to check and 20 understand impacts of the model updates on the simulated mean states. Improvements in the ice water content, the high cloud amounts, the surface air temperature over the Arctic region, the location and the strength of zonal mean subtropical jet, and shortwave radiation over the Africa and the South Asia were found in the NICAM16-S simulations. Some long-standing biases such as the double intertropical convergence zone and smaller low cloud amounts still exist or even worsen in some cases, suggesting further necessity for understanding their mechanisms and upgrading schemes and/or their parameter 25 settings as well as for enhancing horizontal and vertical resolutions.


Introduction
Moist processes play a crucial role in the formation of the Earth's climate. The moist convection redistributes the mass, energy, and momentum of the atmosphere to form large-scale circulation. Clouds are coupled with large-scale circulation evaluations of the mean states simulated by NICAM16-S. Section 5 reports on the computational aspects of the simulation.
Section 6 provides a quick summary of this paper.

Spatial and temporal resolutions 5
Three sets of model configurations were prepared for the HighResMIP simulations and initial and boundary conditions were made for each. NICAM16-S with specific horizontal resolutions were formally labeled NICAM16-7S (56 km mesh), NICAM16-8S (28 km mesh), and NICAM16-9S (14 km mesh). The horizontal mesh size is evaluated as a square root of the mean area of each grid cell (Satoh et al., 2014). The number n in NICAM17-nS is a grid division level (glevel), which denotes a number of subdivisions of the icosahedron to generate a mesh (Tomita et al., 2001). The physics schemes, 10 including parameters, and the initial and boundary conditions are common among different horizontal resolutions except for those explicitly noted in Sections 2 and 3.
The number of vertical levels is 38, with a model top height of around 40 km, equivalent to the previous climate simulations . Atmospheric phenomena of interests may be practically well simulated including the tropical 15 cyclones, MJO, and diurnal precipitation cycle, as we have confirmed in the previous study , though the vertical resolution is not sufficient to resolve the cirrus clouds (Seiki et al., 2015b) and atmospheric gravity waves (Watanabe et al., 2015).
The time step intervals for the dynamical process are set to 240, 120 and 60 s for NICAM16-7S, NICAM16-8S, and 20 NICAM16-9S, respectively. Diffusion coefficients of the divergence damping, the second-order Laplacian background horizontal diffusion, and the first-order Laplacian horizontal diffusion for sponge layer above 20 km are reduced appropriately as the horizontal resolution is increased (Satoh et al., 2008). The time step interval of the cloud microphysics scheme is 30 s and that of the surface processes including turbulence is 60 s for all the horizontal resolutions. The radiation scheme, which requires considerable computational time, is executed every 40, 20, and 10 min for NICAM16-7S, 25 NICAM16-8S, and NICAM16-9S, respectively. Table 1 shows the integration periods of the HighResMIP simulations. For the Tier 1 simulations, which start from 1 January 1950, the initial condition of the atmosphere was taken from the ERA-20C reanalysis (Poli et al., 2016). Strictly following 30 the HighResMIP protocol, the simulations continued until 31 December 2014, using NICAM16-7S and NICAM-8S. The https://doi.org /10.5194/gmd-2019-369 Preprint. Discussion started: 1 April 2020 c Author(s) 2020. CC BY 4.0 License.

Initial conditions
HighResMIP Tier 3 simulations using NICAM16-7S and -8S started from 1 January 2015 as a continuation of the Tier 1 simulations and ended on 31 December 2050. Higher computational cost hinders us from running NICAM16-9S for a hundred years, and thus a time-slice approach is adopted, instead. Specifically, climate simulations have been performed in the following timeframes : 1950-1960, 2000-2010, and 2040-2050. The initial land condition in the past and present-day simulations was taken from a monthly mean climatology of the simulation by NICAM with a mesh size of 220 km (glevel-5) 5 under a present-day condition. It was performed for 10 years, and the last 5 year data were used to obtain the monthly mean land climatology. This approach is the same as the one used in the previous climate simulations ( . The initial land condition for the future time slice run with a 14 km mesh was obtained by interpolating the output of the 28 km mesh run. 10 In addition to the formal HighResMIP simulations, we performed a series of short-term sensitivity experiments to evaluate impacts of the model updates on the simulated climatology, as listed in Table 2. Here, the model configuration of the experimental ID "g" is equivalent to that used in the HighResMIP Tier 1 and 3 simulations, and all the other configurations of the sensitivity experiments are the derivative of "g", as described in Table 2. These experiments were performed from the initial condition beginning from 1 June 2004 because it was used frequently in the previous NICAM studies (e.g., Kodama et 15 al., 2012;Noda et al., 2016;Seiki et al., 2015a). The integration periods of these sensitivity experiments are 1 year in most cases.

External forcings
External forcings of the simulations basically followed the HighResMIP protocol (Haarsma et al., 2016); in other words, 20 historical and SSP585-scenario forcings (O'Neill et al., 2016) were used in the Tiers 1 and 3 simulations.
Daily quarter-degree sea surface temperature (SST) and sea ice mass (ICE) prescribed for the model were obtained from . Since HadISST provides historical sea ice concentration (SIC), ICE was diagnosed from SIC for NICAM (see Section 3.6). Future SIC is estimated from the future SST data and an observed relationship between SST and SIC (Kennedy et al., 2019; https://github.com/PRIMAVERA-H2020/HighResMIP-futureSSTSeaice). Both the SST and ICE were fixed to the boundary conditions in the HighResMIP Tiers 1 and 3 simulations. In the short-term sensitivity experiments, a slab ocean 30 model with a nudging toward the prescribed SST was also tested for a practical purpose (see Section 3.6 for its impact on the simulated climate). In the sensitivity experiments, the fixed SST condition was used in 56 km mesh run whereas and the slab ocean model with the nudging was used in 14 km mesh runs unless explicitly specified. https://doi.org /10.5194/gmd-2019-369 Preprint. Discussion started: 1 April 2020 c Author(s) 2020. CC BY 4.0 License. Figure 1 and Figure 2 display the decadal mean SST and ICE prescribed for the model; Figure 3 exhibits their global-mean variability. Greater warming over the maritime continent, the Indian Ocean, and the edge of the polar regions are found in the 2000s compared with the 1950s, whereas cooling is noticed in the North Pacific and the North Atlantic. Future change in SST is somewhat similar to the El Niño pattern, and the warming is also prominent in the midlatitudes, the tropical Atlantic 5 Ocean, and the edge of the polar regions. Similar tendencies are observed for the distribution of SST trend (not shown). The global mean SST is 17.9 ℃ in the 1950s, 18.2 ℃ in the 2000s, and 19.0 ℃ in the 2040s. ICE continues to decrease from the past to the future. The global mean ICE is anearly halved by the 2040s compared with that in the 1950s.

15
Natural aerosol mass and the number concentration prescribed for the model were obtained from a low-resolution NICAM simulation with online aerosol module based on the Spectral Radiation-Transport Model for Aerosol Species (SPRINTARS) (Takemura et al., 2000(Takemura et al., , 2002(Takemura et al., , 2005(Takemura et al., , 2009Goto et al., 2008Goto et al., , 2011. The simulated climatology of the aerosol in NICAM has been validated (Goto et al., 2018;Suzuki et al., 2008). For CMIP6, NICAM with 220 km mesh (glevel-5) was performed for 100 years using natural aerosol emission with the anthropogenic aerosol module MACv2-SP (Fiedler et al., 2018;Stevens et 20 al., 2017;see Section 3.3) under a perpetual 2012 condition, and the data of the last 90 years were averaged to obtain a monthly-mean climatology of aerosol mass concentration and number concentration of cloud condensation nuclei (CCN) from a natural origin. A low-bound limiter of 50 cm -3 was applied to the CCN prescribed for the model to avoid numerical instability in the cloud microphysics scheme. Figure 4 shows the annual means of the natural aerosol optical thickness and CCN simulated by and prescribed to the model 1 . The climatology of natural aerosol mass and CCN is invariant year by year, 25 whereas anthropogenic aerosol from MACv2-SP is time-dependent in the historical and future simulations (Fiedler et al., 2019;Stevens et al., 2017). Further, aerosol optical properties were overwritten with the stratosphere aerosol dataset (Thomason et al., to be submitted) above the tropopause to introduce effect of volcanic eruptions on the radiation field in a consistent way among different models participating in CMIP6.
As briefly noted in Satoh et al. (2014), a smoother is applied to the model topography to avoid a model abort due to numerical instability. Specifically, a hyper-diffusion is repeatedly applied to the GTOPO30 (doi:10.5066/F7DF6PQS), a global digital elevation model with a horizontal spacing of approximately 1 km, to meet a certain criterion of maximum 10 elevation gradient. The maximum elevation gradient is set to 0.01, 0.01414, and 0.02 mm -1 for NICAM16-7S, NICAM16-8S, and NICAM16-9S, respectively, and is called "A-topography". Note that, in previous NICAM studies using 14 km horizontal mesh (e.g., , "B-topography", in which A-topography of 28 km mesh was interpolated to 14 km mesh grid point, was often used for the sake of stable integration. For CMIP6, A-topography was used to better represent steeper mountains and their effects on the atmospheric phenomena. 15 Table 3 shows a summary of the physics schemes used in NICAM16-S and NICAM.12. Dynamical core and horizontal and vertical diffusion, including the planetary boundary layer schemes in NICAM16-S, are the same as those of NICAM.12 except for some minor updates. A cloud microphysics scheme is used instead of a combination of convection and large-scale 20 condensation schemes to explicitly represent interactions between clouds and circulation. Although most climate models continue to use convection and large-scale condensation schemes even for a mesh size around 14 km, we believe that being free from development and tuning of such complex parameterizations can be another wise option to focus more on investigating the nature of the simulation with explicit cloud microphysics when we expect to proceed toward global cloud resolving climate simulation in the future. Furthermore, global mean precipitation is constrained by radiative cooling in a 25 large-scale clear-sky regions, which can be captured by the relatively coarse resolution run. Although our choice leads to a patchy behaviour of precipitation in the simulation, the simulated climatology of the precipitation pattern, even by the lowest resolution setting (NICAM16-7S), is comparable with the observation (see Section 4).

Model description and impact of model updates on the simulated fields
A series of short-term sensitivity experiments were performed to monitor impacts of several model updates on the simulated 30 climatology. Table 2 offers a full list of sensitivity experiments, and Table 4 presents observations for comparisons with the model . Some of their impacts on the global means are summarized in Table 5, and comparisons with observations along https://doi.org/10.5194/gmd-2019-369 Preprint. Discussion started: 1 April 2020 c Author(s) 2020. CC BY 4.0 License.
with impacts of the horizontal resolutions are shown in Table 6 for reference. Unless explicitly specified, the simulation data are re-gridded to the same grid of observations Table 4 or to 2.5° in latitude and longitude. As noted in Section 3.6, we often prefer to use a slab ocean model with nudging toward the boundary SST rather than the fixed SST condition requested by the HighResMIP protocol because of better performance in the simulated precipitation pattern . Therefore, both the fixed SST and slab ocean configurations were tested in the sensitivity experiments (see Section 3.6). 5

Cloud microphysics
Climate simulations with NICAM16-S were performed with a single moment bulk cloud microphysics scheme with six water categories (hereafter referred to as NSW6). Recently, Roh and Satoh (2014) and Roh et al. (2017) have significantly revised the NSW6 scheme based on comparisons with the TRMM observation. We used the revised version of the NSW6 to 10 show the improvements in the climatology of the NICAM simulation.
The NSW6 scheme was originated with Lin et al. (1983) and Rutledge and Hobbs (1984). After their works, Tomita (2008) modified their cloud microphysics scheme to ensure consistency with thermodynamics used in NICAM; Tomita also simplified it to reduce calculation cost for global high-resolution simulations. The NSW6 was evaluated by comparing the 15 simulated optical properties with satellite observations (Satoh et al., 2010;Kodama et al., 2012;Hashino et al., 2013Hashino et al., , 2016Satoh, 2014, 2018;Roh et al., 2017) using satellite simulators (CFMIP Observational Simulator Package developed by Haynes et al., 2007;Chepfer et al., 2008;Bodas-Salcedo et al., 2008, or Joint Simulator developed by Matsui et al., 2009Masunaga et al., 2010;Hashino et al., 2013). It was revised in each stage of the version management of NICAM (Kodama et al., 2012;Roh and Satoh, 2014;Roh et al., 2017). The revision of the NSW6 scheme by Satoh (2014) 20 andRoh et al. (2017) is the key in this updated version of the NICAM (NICAM16-S). Table 7 summarizes key changes in the NSW6 scheme by Roh and Satoh (2014) and Roh et al. (2017) . In short, the revision aimed to enhance organizations of tropical convective cloud systems by assuming lighter precipitation of graupel and snow and by moderately developing cloud ice. Finally, they achieved successful reproduction of the vertical structures of shallow, 25 congestus, and deep convective clouds over the tropics, comparing to TRMM and CloudSat satellite observations. They also used microwave satellite observations to evaluate the simulated results (Roh and Satoh, 2018); hence, improvements in tropical cloud systems with the revised scheme are robust in terms of optical signals (see the original papers for more details).
Note that a separation of convective and stratiform systems in Roh and Satoh (2014) was omitted in this study because of its small impact (not shown) despite the high computational cost. 30 Table 5a shows the global mean impact of the updated cloud microphysics scheme on the simulated mean values. The most noticeable impact is on an increase in ice water content by more than twice, and this mostly accounts for a snow category in https://doi.org/10.5194/gmd-2019-369 Preprint. Discussion started: 1 April 2020 c Author(s) 2020. CC BY 4.0 License. the cloud microphysics scheme. Figure 5 shows the meridional-height cross section of the observed and simulated zonal mean ice water content (IWC) and its breakdown into the categories of cloud ice, snow and graupel in the simulation. The simulated IWC is largely underestimated using the model before the update of cloud microphysics scheme, as also shown in Seiki et al. (2015a), and it becomes comparable with the CloudSat observation after the update. A noticeable increase in the snow category is seen in the tropical upper troposphere and storm-track region. Cloud ice and graupel are also increased in 5 the upper troposphere, and graupel is rather decreased in the storm-track region. As a result, global mean column-integrated cloud ice and snow are increased by 24 % and 399 % and that of graupel is decreased by 8.4 %. The increase in IWC is consistent with the decelerated development of ice water content by the modified mass and diameter relationship of snow, that reduces the snow density (Table 7d), by the diminished efficiency of accretion of cloud ice by snow (Table 7g), and by ignoring an accretion of snow and cloud ice by graupel (Table 7f). Differences in cloud processes (convection vs. synoptic 10 system) may vary increasing/decreasing in graupel by the model update.
Despite the drastic increase in mass concentrations of snow and cloud ice, the amount of high cloud, particularly, optically thin cloud, is somewhat reduced. Consistently, global mean OLR is increased by about 4 W m -2 (Table 5a), opposite to that found in Roh et al. (2017). In addition, decrease in TOA brightness temperature by the update (Figure 6) is very small 15 compared with that in Roh et al. (2017). These differences between Roh et al. (2017) and this study can be mostly explained by the new coupling procedure between cloud microphysics and radiative transfer as described in Section 3.2 and be partially explained by the different treatment of terminal velocity of cloud ice. The cloud ice terminal velocity diagnosed by Heymsfield and Donner (1990) was replaced with zero in Roh et al. (2017) whereas it was zero and unchanged in this study.
The reduction of the cloud ice fall speed, as in Roh et al. (2017), could increase and elevate the high cloud and decreases the 20 OLR, as also seen in Kodama et al. (2012). A low cloud amount is increased as a result of compensation between an increase in a medium and a thick cloud and a decrease in a thin cloud. These results indicate that the clouds grow thicker on average by updating the cloud microphysics scheme in this study.

Coupling between cloud microphysics and radiative transfer 25
In NICAM16-S, cloud microphysics schemes are fully coupled with a broadband radiative transfer model named MstrnX (Sekiguchi and Nakajima, 2008). The effective radii r e of hydrometeors are calculated with the same assumption of the particle size distribution function as cloud microphysics. including indirect effect, and then passed to the MstrnX. In contrast, fixed effective radii were assumed in NICAM.12 (8 and 40 μm were assumed for liquid and solid hydrometeors, respectively). The following section summarizes the changes in the radiation budget by the coupling between cloud microphysics and radiative transfer as well as details of the update.
The MstrnX requires the database of single scattering properties of hydrometeors (RADPARA), such as the volume 5 extinction coefficient, absorption coefficient, asymmetry factor, and the second moment of phase function (Nakajima et al., 2000). In NICAM.16 we use the RADPARA database revised by Seiki et al., (2014). The RADPARA database of liquid hydrometeors was pre-calculated according to the Mie theory. The non-spherical RADPARA database developed by Fu (1996) and Fu et al. (1998) was applied to solid hydrometeors. The RADPARA database was then compiled as a lookup table of the effective radii from 1 μm to 1 mm to cover size range of most of the hydrometeors in global simulations (Seiki et 10 al., 2014). The effects of precipitating hydrometeors on the radiation budget are detectable specifically over the intertropical convergence zone (ITCZ) and storm-track region (e.g., Chen et al., Li et al., 2014Li et al., , 2016Michibata et al., 2019;Waliser et al., 2011). The revised RADPARA database was evaluated in depth by comparing it with balloon-borne sonde observations in a midlatitude cirrus case (Seiki et al., 2014), and its effectiveness for global simulations was evaluated in several studies (Satoh et al., 2018;Seiki et al., 2015aSeiki et al., , 2015b. 15 Because of non-sphericity, the effective radius of ice particles has a controversial definition, whereas the effective radius is well defined in the case of spherical particles. According to Fu (1996), the effective radius of solid hydrometeors is defined as follows: where = , , are cloud ice, snow, and graupel, is the air density, is the projected area of a particle to flow, and 20 = 916.7 kg m -3 . The integral in equation (1) is analytically calculated using the assumed particle size distribution functions and sponge-like spherical shape in the case of cloud ice and graupel. In contrast, snow has two-dimensional fractal shapes; hence, the numerator-denominator ratio becomes almost constant. Thus, the effective radius of snow is assumed to be constant ( , = 125 μm with = 0.45 2.0 0 ), which is derived by approximating the A-D relationship of aggregates compiled by Mitchell (1996). 25 The impacts of the coupling procedure and non-spherical scattering were examined by comparing a set of seasonal simulations with NICAM-16S (g run in Table 2), NICAM-16S except for the spherical RADPARA database (g9 run in Table 2), and NICAM-16S except for the fixed effective radii and the spherical RADPARA database (g9a in Table 2). Figure 7 shows the zonal mean values of the OLR and the reflected shortwave radiation (OSR) at TOA from the sensitivity 30 experiments. Given the substantial increase in cloud ice and snow from the new NSW6 (cf. Section 3.1), ice optical thickness increases proportionately to the increases in the ice water path (IWP) with the fixed effective radii. As a result, the radiation budget for both the longwave and shortwave is strongly biased in the g9a run. A major portion of the biases is drastically offset in the g9 run by assuming larger effective radii in the coupling procedure. Thus, the coupling procedure automatically prevents artificial biases originating from the inconsistent parameter settings between the cloud microphysics and radiative transfer at the model update. The use of the non-spherical RADPARA database slightly increases OSR over the tropics to the midlatitude of the summer hemisphere because the assumed asymmetry factor for non-spherical particles is 5 smaller than the one for spherical particles (cf. Seiki et al., 2014).
Finally, NICAM16-S shows strong negative biases in OLR over the tropical to subtropical regions and OSR over the polar region and subtropical high-pressure belt. The former biases can be solved by increasing vertical resolution to 400 m near the tropopause with 74 vertical layers (Seiki et al., 2015b). The latter biases mainly stem from the underestimation of low-level 10 clouds since the current updates in cloud microphysics do not work on improvements in warm clouds. The coupling procedure strengthens the negative biases become stronger in the midlatitudes (30°N-60°N) (see g9 and g9a runs in Figure   7) because, in NICAM-12, high clouds associated with extratropical cyclones are artificially brightened and, therefore, conceal the biases due to low-level clouds (cf.

Aerosols in the cloud microphysics and radiation schemes
In NICAM.12, the direct radiative effect of the aerosol is not considered in the radiation scheme, and the number concentration of CCN is set to a constant value of 50 cm -3 , a typical value over the ocean, in the cloud microphysics scheme.
In NICAM16-S, both aerosol's direct and indirect effects are considered by prescribing a distribution of aerosol mass concentration in the radiation scheme and CCN in the cloud microphysics scheme. The natural and stratospheric aerosol data 20 has been described in Section 2.3. For anthropogenic aerosol, a simple plume model, MACv2-SP (Fiedler et al., 2019;Stevens et al., 2017) is used to diagnose the following: (a) a vertical profile of the aerosol optical depth, the single scattering albedo (SSA), and the asymmetry factor and (b) a factor of CCN increase arising from anthropogenic aerosol. This means that the magnitude of the anthropogenic increase in CCN depends on that from the natural origin. Above the tropopause, where volcanic eruptions are major sources of aerosol, the extinction coefficient, SSA, and the asymmetric factor are directly 25 prescribed as CMIP6's external conditions (Section 2.3). Table 5b show an impact of the above updates concerning natural and anthropogenic aerosols on the simulated radiation field. The update reduces the downward shortwave radiation at the surface over most of the continent, particularly over Africa and South Asia, leading to a reduction of an excess of the insolation there (not shown). The reduced insolation at 30 the surface is partly cancelled out by a reduction in the upward longwave radiation over the Africa in association with a decrease in surface air temperature. An enhancement of the surface downward shortwave radiation is dominant over the ocean in association with a thinning of cloud optical depth and a decrease in cloud amount, particularly in the NICAM16-9S simulation (not shown). As a result of these compensations, the global mean net surface radiation change arising from aerosol forcing is around -2 W m -2 in NICAM16-7S and +0.4 W m -2 in NICAM16-9S in this study. Longer integration is needed to assess these impacts quantitatively.

Surface albedo
Surface albedo values were revised based on the observations. In the past, they were tuned to reduce TOA radiation imbalance, which caused a higher bias of surface albedo over the Arctic compared with a satellite observation (Hashino et al., 20 2016). Table 8 shows the revised and previous values of surface albedo used in NICAM. In most cases, the albedo values are set to be smaller in NICAM16-S than those in NICAM.12. In addition to the changes in Table 8, the elevated ocean surface albedo for the direct visible wave by 1.35 times for the radiation calculation in the previous simulations is discarded in NICAM16-S because this is a highly artificial factor.

25
Consistent with the reduced albedo and the increased net upward shortwave radiation, sensitivity experiments using NICAM16-7S show that the use of the new albedo configuration tends to reduce the surface air temperature bias over the land ice compared with the old one (Figure 10a vs. Figure 10b and black vs. green lines in Figure 10d). Specifically, the cold bias in Greenland, the Himalayas, and the Antarctic is reduced. Upward net longwave radiation at the surface is decreased over the ocean, consistent with the increased surface albedo for infrared (i.e., less blackbody) over the ocean. Warm bias still 30 exists in the Arctic in the new configuration, and this will be reduced by changing a sea ice configuration, as explained in Section 3.6.

Treatment of ocean
A mixed-layer slab ocean model similar to McFarlane et al. (1992) had been implemented in NICAM. The model predicts SST, ICE, snow over sea ice, and snow temperature by solving a heat balance between ocean, sea ice, snow, and the atmosphere. A depth of the slab ocean model is set to 15 m, considering the better performance of the simulated precipitation 5 pattern  and MJO (Grabowski, 2006). A simple nudging technique is used to force the predicted SST toward a reference SST with a relaxation time of 7 days. In the slab ocean model implemented in NICAM.12 and NICAM.16, SIC is diagnosed from ICE, as, where SICCRT is set to 300 kg m -2 .

10
In many cases including the HighResMIP protocol, only the SST and SIC data are provided from CMIP6 to run the model, and, therefore, the ICE data prescribed for the model should be diagnosed from SIC. In NICAM.12 and NICAM16-S, ICE is diagnosed simply as, In the previous study using NICAM.12, we often set a value of SICCRT to 300 kg m -2 for Eq. (3), the same as that used in Eq. (2). However, this situation causes an underestimation of ICE over most of the sea ice areas and leads to a warm bias 15 over the Arctic . Based on an ocean model result (Tatebe, personal communication) and sensitivity experiments (Figure 10c), SICCRT is set to 1,600 kg m -2 for Eq. (3) in NICAM16-S to diagnose ICE from SIC.
In the HighResMIP protocol, SST and SIC in the model are fixed to the time-varying external conditions. This fixed SST/SIC condition is achieved by nudging SST and ICE toward the prescribed external condition with a zero-relaxation time 20 in the slab ocean model. However, fixed SST simulation is known to cause severe bias in the precipitation pattern in the tropics , and thus the use of the slab ocean model with a 7 day relaxation time is often preferred. Table   5e and Figure 11 summarize the comparisons of the simulations between the fixed SST condition and the slab ocean model and nudging. Overall, the global mean impact of the slab ocean model is not very large. In the simulation with slab ocean model, global mean precipitation shows a slight increase and OLR shows an increase, which is associated with a slightly 25 warmer surface air temperature. The introduction of the slab ocean model considerably affects the horizontal distribution of the cloud and precipitation system. Double ITCZ bias is more prominent in the precipitation, high cloud fraction, and OLR fields in the fixed SST runs compared with the slab ocean runs, particularly in the high-resolution run. As far as our investigation shows, NICAM16-9S with the slab ocean model best simulates the ITCZ peak precipitation and precipitation pattern. Although an importance of the short-term SST variation driven by the atmosphere on the pattern of the precipitation is apparent from Figure 11, the introduction of the slab ocean model alone does not resolve the bias of the precipitation pattern. Further analysis is necessary to understand the physical mechanisms of the bias, notably perhaps the timescale of the convection,. 5

Orographic gravity wave drag
In NICAM16-S, the conventional orographic gravity wave drag scheme (McFarlane, 1987) is tested to better simulate the location and strength of the subtropical jet. The wave generation parameter (α), which is proportional to the product of wave generation efficiency and representative horizontal wavenumber (Eq. 3.1b in McFarlane 1987), is roughly doubled as the 10 horizontal mesh size is halved. Specifically, α is set to 3.38 × 10 −5 for 56 km mesh, 7.12 × 10 −5 for 28 km mesh, and 1.46 × 10 −4 for 14 km mesh, respectively. Figure 12 and Figure 13 show the zonal mean zonal wind in boreal summer and winter, simulated with and without the gravity wave drag scheme. The zonal mean zonal wind simulated without the gravity wave drag scheme is biased poleward in both 56 km and 14 km mesh runs. The gravity wave drag scheme decelerates the zonal mean zonal wind at the poleward flank of the subtropical jet, especially in NH winter, reducing the locational bias of 15 the jet. The impact of the gravity wave drag is larger in the 56 km mesh model than that in the 14 km mesh model. The pattern of the response of the zonal mean zonal wind to the orographic gravity wave drag scheme is similar to that of previous studies (e.g., Iwasaki et al., 1989;McFarlane, 1987).
Although it is believed that even a 14 km mesh is insufficient to explicitly simulate the effects of the orographic gravity 20 wave drag on the mean field (Nappo, 2012), it may not be a wise choice to introduce such a gravity wave drag scheme to the global non-hydrostatic model. Gravity wave drag scheme introduces the uncertain parameter α, which is tuned to best simulate the climatology of the zonal wind for each resolution, although we did not tune α for each resolution. There is no solid guideline to determine α, including the dependency of the wave generation efficiency and representative horizontal wavenumber on the horizontal and vertical resolutions. Therefore, the use of a gravity wave drag scheme may hinder the 25 pure resolution dependency of the mean field and suppress the advantages of the high-resolution model for simulating largescale circulation in a seamless manner. Nevertheless, we determined to use the orographic gravity wave drag scheme for CMIP6 to reduce the locational bias of the subtropical jet so that an improvement of the tropical cyclone track would ensue.
It is important to recognize the merits and demerits involved in the use of a gravity wave drag scheme and reconsider its use depending on the main purpose of the simulation. Table 6 shows the summary of the comparisons between NICAM16-S simulations and the observations along with a dependency of horizontal resolution. In terms of the global mean energy budget of the atmosphere, most of the bias features are similar to the previous NICAM climate simulations with 14 km mesh , meaning that both OLR and OSR are underestimated. In the past, a high cloud amount was overestimated in the NICAM simulations. Now, a high cloud 5 amount is comparable with the ISCCP observation in terms of global mean, although higher altitudes and thinner optical depth (not shown) are simulated, leading to the underestimation of OLR. Global mean precipitation and OLR are decreased as the resolution is increased, consistent with a previous study using 3.5-14 km mesh NICAM (Miyakawa and Miura, 2019).

Preliminary evaluations with observations including dependency of horizontal resolution
Low-level cloud amount is substantially underestimated, especially in the NICAM-7S and NICAM-8S runs, leading to the underestimation of OSR. Surface air temperature is slightly decreased as the resolution is increased in association with an 10 increase in total cloud fraction. The most noticeable change from the previous simulation in terms of the global mean is the IWP. As described in Section 3.1, IWP is drastically increased by the update of cloud microphysics scheme. Table 9 shows the elapsed time and the file size of the simulations by NICAM16-7S, NICAM16-8S and NICAM16-9S on the Earth Simulator 3 (NEC SX-ACE). File staging option was used in the NICAM16-8S and NICAM16-9S simulations, whereas the file was directly read from the global file system in the NICAM16-7S simulation to save time for queuing. The actual wall clock time, including queuing time, was a few times greater than SYPD in Table 9 for NICAM16-8S and NICAM16-9S depending on the congestion of the Earth Simulator 3. For the product run of NICAM16-8S, 160 nodes are 20 used to finish a 101 year simulation within a realistic time.  (Yashiro et al., 2016), the measurement included an initial setup and input/output processes. Physics is a major part that contributes to the total elapsed time. Among the several physics components, the radiation scheme primarily contributes 25 to the total elapsed time, followed by the cloud microphysics scheme, consistent with Yashiro et al. (2016) on the K computer. As the resolution increases, the percentage of the dynamics is increased and that of the cloud microphysics and near-surface processes decreases because time interval of the cloud microphysics and near-surface processes is set to be invariant throughout the models with different resolutions. Because the dynamics involves communication at every time step, some parts of the elapsed time counted as the dynamics can actually be waiting time because of node imbalance occurring in 30 other processes.

Post-processes
The data are output in an icosahedral grid on the height above sea level (ASL) or on the standard pressure. The vertical interpolation from terrain-following height to the ASL height or standard pressure is performed online using second-order Lagrange interpolation during the simulation. 5 Post-processes are performed in the following order: ico2ll, roughen, and z2pre.
(1) ico2ll All the native icosahedral grid data are converted to high-resolution latitude-longitude grid data by area-weight averaging.
The interval of latitude and longitude is determined so that the longitudinal interval is close to the average interval of the 10 icosahedral grid (Satoh et al., 2014) on the equator. Specifically, the interval of longitude and latitude is 0.56° for NICAM16-7S, 0.28° for NICAM16-8S, and 0.14° for NICAM16-9S.
(2) roughen The high-resolution latitude-longitude grid data are coarsened to low-resolution latitude-longitude grid data by area-weight averaging. It is often necessary to reduce the data size by coarsening, although CMIP6 does not request us to coarsen the 15 data. We prepare 1.0°, 1.25°, and 2.5° data for analysis.
(3) z2pre Several three-dimensional variables are converted from the ASL height to the standard pressure at this point by linear interpolation, if it is necessary.
Pressure velocity is diagnosed from the vertical velocity and temperature after ico2ll, assuming hydrostatic balance.
Geopotential height is calculated from the linear interpolation of vertical levels using logarithms of pressure, assuming constant gravity acceleration irrespective of height, which is a consistent treatment with the model configuration. 25

Summary
This paper described the experimental design, model description and impacts of the model updates on simulated climatology using NICAM prepared for CMIP6 HighResMIP (NICAM16-S). The major updates and their impacts are summarized as follows: • Update of the cloud microphysics scheme: Snow and cloud ice increase in the atmosphere, leading to less high cloud amount and more OLR.
• Implementation of the coupling between cloud microphysics and radiation schemes: The negative OLR bias is reduced in association with the larger cloud ice effective radius.
• Update of treatment of natural and anthropogenic aerosols: Local surface radiation budget is improved, especially over 5 Africa and South Asia.
• Update of land model: Overall, the soil moisture and the precipitation increase and the surface air temperature decreases over the continent.
• Revision of the surface albedo values: Cold bias in the Greenland, the Himalayas and the Antarctic is reduced.
• Change in the diagnostics of ICE: Warm bias over the Arctic region is reduced. 10 • Introduction of the orographic gravity wave drag scheme: The location and strength of the zonal mean jet are improved.
Comprehensive evaluations and future projection using full HighResMIP data by NICAM16-S will be presented in a forthcoming paper.

Code and data availability 15
The exact model source code, input data or scripts to generate them, and scripts for the simulations and the post-processes used to produce the results presented in this paper are archived on Zenodo (doi:10.5281/zenodo.3727329). The model source codes are shared by the NICAM community and available for those who are interested as long as a user follows the terms and conditions described in http://www.nicam.jp/hiki/?Research+Collaborations. Most of the input data are freely accessible from input4MIPs (https://esgf-node.llnl.gov/projects/input4mips/) for ocean boundary condition, GHG concentration, ozone 20 and solar forcing, from ECMWF website (https://apps.ecmwf.int/datasets/data/era20c-daily/) for ERA-20C reanalysis, from

Author contributions
CK and ATN managed overall CMIP6 activity in NICAM group and prepared the initial and boundary conditions, and MS managed development and scientific activity in NICAM group. CK added interfaces of the initial and boundary conditions to NICAM, and TO added a function to output variables requested by CMIP6 and converted the output data using CMOR3. CK, TO, TS, HY, MN, YY contributed to the development of NICAM16-S including bug fixes, and WS, TN, DG provided their 5 schemes and parameters for the development. CK performed the product runs and all the sensitivity experiments described here, transferred the data to ESGF, and wrote a major part of this paper. TS wrote most part of Section 3.1 and 3.2 and TS, ATN, DG, HM, TN modified the manuscript. All the authors provided advice for the development of NICAM16-S and/or experimental design and reviewed the manuscript.

Competing interests
The authors declare that they have no conflict of interest.

Acknowledgment
The authors would like to thank Hiroaki Tatebe and Ryosuke Shibuya for discussion on model configurations and Manabu 15 Abe and Takahiro Inoue for technical advice for CMIP6. CERES, CloudSat, ISCCP, and GPCP data are obtained from the National Aeronautics and Space Administration (NASA), JRA-55 data from the JMA, and GridSat data from the National Goto, D., Nakajima, T., Takemura, T. and Sudo, K.: A study of uncertainties in the sulfate distribution and its radiative forcing associated with sulfur chemistry in a global aerosol model, Atmos. Chem. Phys., 11(21), 10889-10910, 20 doi:10.5194/acp-11-10889-2011

Descriptions g
Same as NICAM16-S.

g3
Same as g but for using the previous cloud microphysics scheme used in NICAM.12 (Table 7; Section 3.1).

g6
Same as g but for prescribing zero natural and anthropogenic aerosol mass concentration for the radiation scheme and constant CCN of 50 cm -3 for the cloud microphysics scheme (Section 3.3).

g4
Same as g but for omitting the effects of wetland and accumulation of water on land ice (Section 3.4).

g8
Same as g but for switching off the subgrid-scale orographic gravity wave drag scheme (Section 3.7).

g9
Same as g but for considering only the spherical particle in the radiation table (Section 3.2).

g9a
Same as g9 but for removing the interaction between radiation and cloud microphysics (Section 3.2).

f1d
Same as g but for a previous configuration of aerosol effective radii # .

f1
Same as f1d but for using the previous SICCRT value of 300 kg m -2 (Section 3.6).
f Same ad f1 but for using the previous surface albedo values (Table 8; Section 3.5).
# In the f1d run, the aerosol effective radii of soil dust and seasalt are set to 4 μm and 2 μm, respectively, following the older MstrnX (Nakajima et al., 2000). In the g (=NICAM16-S) runs, those of soil dust and seasalt are set to 1.6 μm (Mahowald et al., 2014;Omar, 2005) and a function of relative humidity (considering hygroscopity).       following Thompson et al., (2008).
Ice hydrometeors are assumed as the spherical shape with fixed bulk densities following Rutledge and Hobbs (1984). # The particle size distribution of rain was also revised in the original paper, but the revision is not used in the latest version because of computational efficiency against its small impact on improvements. In addition, the assumption that cloud ice does not precipitate is inconsistent with some ice cloud microphysics but is used to reproduce observed high cloud signals 5 over the tropics ad hoc.       Table 7; g3 and g runs in Table 2). Breakdown of the simulated IWC into cloud ice (d), snow (e), and graupel (f) is shown on the bottom panels. The contour shows the simulation using the cloud microphysics schemes before the update and shading shows the difference between 5 the simulations using the cloud microphysics schemes before and after the update. The analysis data are 0.25° (a) and 2.5° (b-f) gridded data, and the vertical axis is the altitude above sea level.