LICOM Model Datasets for the CMIP6 Ocean Model Intercomparison Project

The datasets of two Ocean Model Intercomparison Project (OMIP) simulation experiments from the LASG/IAP Climate Ocean Model, version 3 (LICOM3), forced by two different sets of atmospheric surface data, are described in this paper. The experiment forced by CORE-II (Co-ordinated Ocean–Ice Reference Experiments, Phase II) data (1948–2009) is called OMIP1, and that forced by JRA55-do (surface dataset for driving ocean–sea-ice models based on Japanese 55-year atmospheric reanalysis) data (1958–2018) is called OMIP2. First, the improvement of LICOM from CMIP5 to CMIP6 and the configurations of the two experiments are described. Second, the basic performances of the two experiments are validated using the climatological-mean and interannual time scales from observation. We find that the mean states, interannual variabilities, and long-term linear trends can be reproduced well by the two experiments. The differences between the two datasets are also discussed. Finally, the usage of these data is described. These datasets are helpful toward understanding the origin system bias of the fully coupled model.


Background and summary
The Ocean Model Intercomparison Project (OMIP) is one of the endorsed Model Intercomparison Projects in phase 6 of the Coupled Model Intercomparison Project (CMIP6). The purpose of OMIP is to understand the origin of systematic model biases. The subject also belongs to one of the scientific questions of CMIP6. The primary contribution of OMIP to the World Climate Research Programme's Grand Challenges in Climate Science is in regional sea level change and near-term (climate/decadal) prediction. The OMIP experiment is a hindcast experiment for global coupled ocean-sea-ice models, inheriting the protocol from the Co-ordinated Ocean-Ice Reference Experiments, Phase II (CORE-II; Griffies et al., 2009). The experiments are driven by the prescribed atmospheric forcings from modified reanalysis and observational datasets. The turbulence momentum and heat fluxes are computed through the method of Large and Yeager (2004).
The LASG/IAP Climate Ocean Model (LICOM) was developed at the State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics (LASG), Institute of Atmospheric Physics (IAP), Chinese Academy of Sciences (CAS) around 30 years in the late 1980s (Zhang and Liang, 1989). The latest version of LICOM is LICOM3, preparing for CMIP6. Currently, LICOM3 is the ocean component of two climate system models of CMIP6: the CAS Flexible Global Ocean-Atmosphere-Land System model, finite-volume version 3 (CAS FGOALS-f3) and the CAS Flexible Global Ocean-Atmosphere-Land System model, grid-point version 3 (CAS FGOALS-g3), which in CMIP5 were FGOALS-s2 (Bao et al., 2013, Lin et al., 2013a and FGOALS-g2 (Li et al., 2013;Lin et al., 2013b), respectively. In addition, the previous version of LICOM, LICOM2.0 (Liu et al., 2012), is also employed as the ocean component in CAS-ESM (the CAS Earth System Model), version 1.0 (Dr. Minghua ZHANG, personal communication, 2016).
The OMIP simulations of LICOM3 were finished in June 2019 and the data have submitted to the Earth System Grid (ESG) data server (https://esgf-nodes.llnl.gov/projects/ cmip6/). The purpose of the present paper is to provide a comprehensive description of the OMIP datasets of LICOM3 for a variety of users. We document detailed descriptions of the model configurations and experiments, as well as the algorithm of the diagnostic variables. Section 2 presents the model descriptions and experiment designs. Section 3 presents a basic technical validation of the LICOM3 experiments. Section 4 describes the datasets. Section 5 provides usage notes.

Model and experiments
2.1. Introduction to the model Since LICOM2.0 (Liu et al., 2012), LICOM has been substantially upgraded in the interface with the flux coupler, the dynamic core, and the physical packages (Table 1). First, we upgraded the LICOM interface from the NCAR flux coupler 6 to coupler 7 (Lin et al., 2016), because the new version has been optimized for high-resolution modeling (Craig et al., 2012). Here, LICOM3 coupled with the Community Ice Code, version 4 (CICE4), i.e., the ocean-ice coupled model, is used to conduct the OMIP experiments. The prescribed atmospheric data have been input from the at-mospheric data model and then passed to the coupler to drive the ocean-sea-ice coupled model. Second, the orthogonal curvilinear coordinate (Madec and Imbard, 1996;Murray, 1996) has been introduced in LICOM. The tripolar grid can be used with the North Pole split into two poles on the land in the Northern Hemisphere (NH) at (65°N, 65°E) and (65°N, 115°W). This improvement has solved the problem of the singularity of the North Pole in the normal longitude-latitude grid. Meanwhile, the spatial filter in the high latitudes of LICOM has also been eliminated, and the scalability and efficiency of the parallel algorithm has been extensively improved. Besides, preserved shape advection in the tracer formulation (Xiao, 2006) and implicit vertical viscosity (Yu et al., 2018) are also employed in the present version.
Third, the tidal mixing of St. Laurent et al. (2002) and the buoyancy frequency (N 2 )-related thickness diffusivity of Ferreira et al. (2005) have been introduced into LICOM after CMIP5. The effects of tidal mixing on the Atlantic meridional overturning circulation (AMOC) was preliminary evaluated by Yu et al. (2017). The effects of thickness diffusivity with temporal and spatial variation were evaluated by Li et al. (2019). Besides, the chlorophyll-a-dependent solar penetration of Ohlmann (2003), the vertical diffusivity of Canuto et al. (2001Canuto et al. ( , 2002, and the isopycnal mixing of Redi (1982) and Gent and McWilliams (1990, GM90 hereafter) are also used in LICOM3. The isopycnal mixing coefficient is constant, with a value of 300 m 2 s −1 . The thickness mixing coefficient (i.e., diffusivity) for GM90 employs the scheme of Ferreira et al. (2005), in which the coefficient is dependent on the spatial distribution of N 2 and varies with location and time. The coefficient is set to 300 m 2 s −1 within the mixed layer or in the coastal region where the water depth is shallower than 60 m, while it varies between 300 and 2000 m 2 s −1 in other places.
The virtual salinity flux is computed as the freshwater Mixing in the mixed layer (Canuto et al., 2001(Canuto et al., , 2002 Mixing in the mixed layer (Canuto et al., 2001) and internal tide mixing (St. Laurent et al., 2002) Isopycnal mixing Isopycnal mixing (Redi, 1982) and advection (Gent and McWilliams, 1990) Isopycnal mixing (Redi, 1982) and advection (Gent and McWilliams, 1990) with N 2 thickness diffusivity (Ferreira et al., 2005)  flux multiplied by a constant salinity of 34.7 psu. A restoring term with a piston velocity of 20 m yr −1 has been applied to the virtual salinity flux. If there is sea ice, a piston velocity of 50 m (20 d) −1 is applied under sea ice. Here, the low-resolution LICOM3, which is used both for CMIP6 and OMIP, has 360 and 218 grid numbers for the zonal and the meridional directions, respectively. LICOM3 has two resolutions in the vertical (30 and 80 levels), but only the 30-level resolution is used for OMIP and CMIP6 to save computing resources. The depths of the W-grid and T-grid are shown in Table 2. At the same time, an eddy-resolving version of LICOM3 has also been developed, with a horizontal resolution of about 10 km at the equator and 2.7 km around the Antarctic, and 55 levels in the vertical. This version has also been implemented in an ocean forecast system for short-term ocean prediction 1 .

Experiment designs
The OMIP experiment follows the protocol of CORE-II, which is an ocean-sea-ice coupled hindcast simulation forced by about 60 years of modified reanalysis atmospheric variables with diurnal to decadal signals. Usually, the experiment is conducted for five or six cycles to reach an equilibrium state. Details of the experiments are shown in Table 3. Here, LICOM3 is coupled with CICE4 using the NCAR flux coupler 7. Two standard OMIP experiments have been conducted: one forced with CORE-II data derived from NCEP-NCAR reanalysis (Large and Yeager, 2004), named OMIP1; and the other forced with the surface dataset for driving ocean-sea-ice models based on Japanese 55-year atmospheric reanalysis (JRA55-do, Tsujino et al., 2018), named OMIP2. The atmospheric variables include atmospheric surface wind vectors at 10 m, temperature at 10 m, specific humidity at 10 m, air density at 10 m, precipitation (both rain and snow), surface downward shortwave radiation, downward longwave radiation, and sea level pressure. The period of the forcing data is 62 years  for CORE-II, while the period is 61 years (1958-2018) for JRA55-do. Although the frequencies of the two datasets are different, we use 6-h intervals for both datasets to force the model. Here, experiments with six cycles (one cycle corresponds to 1948-2009 for CORE-II and 1958−2018 for JRA55-do) have been completed and uploaded onto ESG nodes for both OMIP1 and OMIP2. Details of the primary output and diagnostic variables are given in Table 4.

Global mean variables
The basic results of OMIP1 and OMIP2 from LICOM3 are evaluated before submitting the whole datasets. Here, some primary aspects are selected to demonstrate the performance of the model. After six cycles of integration, the sur-face fluxes are balanced by the adjustment of intrinsic processes and the model reaches a quasi-steady state. Figure 1 shows six cycles of the global mean sea surface temperature (SST), the volume mean temperature, the global mean sea surface salinity (SSS), and the volume mean ocean temperature for both OMIP1 (blue curve) and OMIP2 (red curve). The global mean SST value quickly reaches an equilibrium state (Fig. 1a). The global volume mean temperature increases at a rate of 0.09°C (100 yr) −1 due to the input of net heat flux at the sea surface (Fig. 1b). The gradual decrease in net heat flux (from 0.4 W m −2 for the first cycle to 0.17 W m −2 for the sixth cycle), due to the SST increase, leads to the decrease in the temperature trend from cycle to cycle. The simulated global mean temperature reaches a quasi-stable state after the fifth to the sixth cycle. The global mean SSS value has a small trend, with a value of 0.02-0.03 psu (100 yr) −1 , and reaches a near steady state after the second cycle. Meanwhile, the global mean volume ocean salinity has a trend of 0.01 psu (100 yr) −1 due to the positive virtual salt flux at the sea surface. The global annual mean SSTs from the last cycle of the two experiments are compared with Extended Reconstructed SST, version 5 (ERSST.v5) data (Huang et al., 2017), which serve as the observed reference values (Fig. 2a)  Associated with global warming, the sea-ice areas in the polar regions also shrink significantly. The change in sea ice will further affect the albedo and greatly perturb the heat entering the ocean. The simulated annual mean NH and Southern Hemisphere (SH) sea-ice cover (SIC) are shown along with their observed values from the National Snow and Ice Data Center (NSIDC; Fetterer et al., 2017;Fig. 2b). The annual mean values in the NH for models and observations match very well. The correlation coefficients between OMIP1 and OMIP2 and observations during 1980-2009 are 0.89 and 0.96, respectively. The simulated amplitude of interannual-decadal variability for OMIP2 is also better than that for OMIP1. The STDs are both 0.2 × 10 6 km 2 for NSIDC and OMIP2, but much larger for OMIP1 (0.34 × 10 6 km 2 ).
In the SH, the results from OMIP2 are also much closer to observation than that of OMIP1. The annual mean SIC val-ues are overestimated by about 1.0-2.0 × 10 6 km 2 in OMIP1 and OMIP2, although the latter is closer to the observation. The correlation coefficient between the simulation and the observation is 0.81 for OMIP2, but only 0.41 for OMIP1, during 1980-2009. The STDs are 0.24, 0.47 and 0.26 × 10 6 km 2 for NSIDC, OMIP1 and OMIP2, respectively.
The simulated linear trends of SIC are negative for both the NH and SH. However, both OMIP1 and OMIP2 underestimate the magnitude, which is similar to the situation for SST. The observed linear trends are −0.50 and 0.15 × 10 6 km 2 (10 yr) −1 for the NH and SH, respectively. For comparison, the simulated trends are −0.48 and 0.08 × 10 6 km 2 (10 yr) −1 for OMIP1 and −0.38 and 0.14 × 10 6 km 2 (10 yr) −1 for OMIP2. The global mean thermosteric sea level and the volume transport at 26.5°N in the Atlantic, which is used to meas-   (SST during 1958-2009, SIC for 1980-2009, 2005for Argo, 2004 for RAPID) are provided in the figures. ham et al., 2007). The global mean thermosteric sea level reflects the heat uptake of the ocean. Although the simulation can reproduce the interannual and interdecadal variations of the whole system, the increasing/decreasing trend magnitudes in the model are also smaller than those in the observation, which is consistent with the trends of both SST and SIC. The simulations of AMOC for OMIP1 and OMIP2 are larger than that from RAPID, although the simulated values for OMIP2 are closer to the observation. The simulations for both OMIP1 and OMIP2 can capture the weakened AMOC from 2005 to 2009, but the simulation for OMIP2 fails to reproduce the increasing AMOC transport after 2012. This bias can also be found in other OMIP2 simulations, indicating a possible issue of OMIP2 forcing in the northern Atlantic Ocean (Dr. Hiroyuki TSUJINO, personal communication, 2019).

SST and SSS biases
To understand the origin of systematic biases of ocean model, the biases of SST and SSS are shown in Figs. 3a-b and Figs. 3d-e, respectively. In general, the SST biases for the two experiments are similar, including large cold SST biases located in the regions of strong western boundary currents and the Arctic connected with the North Atlantic, and warm SST biases in the eastern boundary regions. The rootmean-square errors (RMSEs) of SST are 0.66°C for OMIP1 and 0.63°C for OMIP2. The improvements in SST simulation for OMIP2 occur in the eastern boundary, the warm pool, southern Indian Ocean, the Norwegian Sea, and Barents Sea (Fig. 3c).
The SSS for OMIP1 and OMIP2 have the same bias patterns and the same RMSEs (0.45 psu). The common biases include larger than 1.5 psu salty biases in the East Siberian Sea, Chukchi Sea and Beaufort Sea of the Arctic Ocean; The improvements in SSS simulation for OMIP2 occur in the South China Sea, eastern Indian Ocean, western Atlantic Ocean, and in the Southern Ocean (south of 40°S; Fig. 3f). However, the SSS biases for OMIP2 increase in the western Indian Ocean. In the eastern Indian Ocean and western Pacific warm pool near the western coast, the salty SSS bias becomes the fresh SSS bias from OMIP1 to OMIP2. Meanwhile, the salty biases for OMIP2 also increase slightly in the subtropics of the South Pacific.

Sea surface height
The simulated sea surface height (SSH) is compared with satellite data from AVISO (Archiving, Validation and Interpretation of Satellite Oceanographic Data; http://www. aviso.altimetry.fr/) in Fig. 4. The observed high SSH is located in the western Pacific and Atlantic, eastern Indian Ocean, and southwest Indian Ocean, while the low SSH is located in the Southern Ocean, Arctic Ocean, and sub-polar gyres in the high-latitude Atlantic and Pacific. The model can capture all the high or low features very well. The spatial correlations between the simulations and observation are about 0.98. This indicates that the upper-layer horizontal circulations are also simulated well.

AMOC
The AMOC results are displayed in Fig. 5. Both OMIP1 and OMIP2 can capture the features of that derived from the World Ocean Circulation Experiment hydrograph-  ic sections (Lumpkin and Speer, 2007). In the upper ocean (0-400 m), there are wind-driven cells connecting the tropics and subtropics. The cell is northward (southward) in the NH (SH). In the NH, the cell includes an equatorial current system and western boundary currents (Gulf Stream). In the middle and high latitudes of the NH, the upper-ocean cell is southward, corresponding to the sub-polar gyre. The maximal volume transport of North Atlantic Deep Water (NADW), which is located at the latitude of 38°N and a depth of about 1000 m, can reach 22 Sv. The upper branch of NADW can enter the Arctic Ocean (> 60°N). The low branch of NADW returns southward. The minimal volume transport of Antarctic Bottom Water (AABW) is −6 Sv and limited below 3500 m.
Compared with the observed one at 26.5°N during the period 2005-09, the depth of maximal AMOCs is close to the observed. The depth from NADW shifting to AABW is located at about 3000 m for both simulations, which is shallower than the 4500 m for the observed and robust diagnostic values from Lee et al. (2019). The upward shift of the transition depth is a common problem for surface forced oceanice simulations (Danabasoglu et al., 2014(Danabasoglu et al., , 2016. The simulated magnitudes of maximal AMOC transport are likely overestimated by 0.5-2.2 Sv than observation at 26.5°N. At 26.5°N, the AMOC for OMIP2 is about 18.24 Sv, which is slightly larger than the observed value from RAPID (17.57 Sv), but smaller than the value for OMIP1 (20.75 Sv).

Interannual-decadal SST amplitude
The interannual-decadal SST amplitude is denoted as the STD. The STDs of SST are presented in Fig. 6. The simulated amplitudes from OMIP1 and OMIP2 can capture well the observed one, with spatial correlations of 0.96. The large amplitudes are located in the tropical Pacific, the North Atlantic (> 40°N) and North Pacific (> 30°N) for the observation and simulations. However, the amplitudes are overestimated in the tropical Pacific and North Atlantic both in OMIP1 and OMIP2. Further comparisons show the simulated amplitude is closer to the observed in OMIP2 than that in OMIP1, including the tropical Pacific and North Atlantic. In the Niño3 region, the amplitudes are overestimated by about 32% and 16% for OMIP1 and OMIP2, respectively, compared with the observed (0.91). This indicates the OMIP2 forcing is better for simulating interannual-decadal variabilities than that of OMIP1.

Data records
The two sets of data, OMIP1 and OMIP2, have been uploaded onto ESG nodes and can be found at https://esgfnodes.llnl.gov/projects/cmip6/. The dataset format is Network Common Data Form (NetCDF), version 4. Although the model outputs are double precision, we converted all the variables into single precision for analysis. These data can be easily handled by common computer programming languages and professional software, such as Climate Data Operators (https://code.mpimet.mpg.de/projects/cdo/) or NetCDF Operator (http://nco.sourceforge.net).

Usage notes
The original model outputs are on a tripolar grid with two poles set on land in the NH. The horizontal grid numbers are 360 and 218 in the zonal and meridional directions, respectively. The original grid distribution is kept and the format is slightly changed to CMOR (Climate Model Output Rewriter) file structure as required by OMIP. The data have 30 vertical levels and the original vertical level is not changed on ESG nodes. The first level is at the depth of 5 m with the thickness of 10 m. The depths of each level can be found in Table 3.
The variables of Priority 1 for OMIP include monthly ocean temperature, salinity, ocean velocities, SSH, etc. (Table 4). There are 31 Priority 1 variables. Additionally, the daily maximal ocean mixed-layer thickness, defined as the depth that potential density is 0.03 kg m −3 larger than that on the sea surface, is also provided. Data grid information, such as area, mask, cell thickness and volume for every ocean grid, is also provided.