K-1 Coupled GCM (MIROC)
Description

K-1 model developers
Center for Climate System Research (CCSR), University ofTokyo; National Institute for Environmental Studies (NIES);
Frontier Research Center for Global Change (FRCGC)

Edited by Hiroyasu Hasumi and Seita Emori

September, 2004

The K-1 project, formally the subject No. 1 of the Kyousei project, is a part of the
Human-Nature-Earth Symbiosis Project. It is supported by the Research Revolution 2002
Program of the Ministry of Education, Culture, Sports, Science, and Technology of Japan.
K-1 model developers (in alphabetical order):ABE-OUCHI, Ayako (CCSR/FRCGC)EMORI, Seita (NIES/FRCGC)
HASEGAWA, Akira (FRCGC)
HASUMI, Hiroyasu (CCSR)INOUE, Takahiro (Research Organization for Information Science and Technology)
KIMOTO, Masahide (CCSR)
MATSUMURA, Shinji (FRCGC)NAGASHIMA, Tatsuya (NIES)
NAKANO, Hideyuki (CCSR, currently at Meteorological Research Institute)
NISHIMURA, Teruyuki (FRCGC)NOZAWA, Toru (NIES)OHGAITO, Rumi (FRCGC)
OKA, Akira (CCSR)OKADA, Naosuke (NIES)
OGOCHI, Koji (NEC Solutions)
OGURA, Tomoo (NIES)
SAITO, Fuyuki (CCSR)SAITO, Kazuyuki (FRCGC)
SAKAMOTO, Takashi (FRCGC)SEGAWA, Tomonori (FRCGC)
SOGA, Takashi (NEC Solutions)
SUDO, Kengo (FRCGC)
SUMI, Akimasa (CCSR)SUZUKI, Tatsuo (FRCGC)SUZUKI, Tsuneaki (FRCGC)
TAKAHASHI, Hideki (NEC Solutions)
TAKATA, Kumiko (FRCGC)
TAKEMURA, Toshihiko (Kyushu University)
TAKIGAWA, Masayuki (FRCGC)
TSUSHIMA, Yoko (FRCGC)WATANABE, Masahiro (Hokkaido University)
WATANABE, Shingo (FRCGC)
YOKOHATA, Tokuta (NIES)

We would like to express our deepest gratitude to late Dr. Atusi Numaguti, who was theprincipal developer of CCSR/NIES AGCM and devised the framework for the ﬂux couplerof our model.
Corresponding to:
Hiroyasu Hasumi
Center for Climate System Research, University of Tokyo
4-6-1 Komaba, Meguro-ku, Tokyo 153-8904, Japan
E-mail: hasumi@ccsr.u-tokyo.ac.jp

Contents
1 Overview: Structure of MIROC

1

1.1

Flux coupler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Coupling strategy for parallel computing . . . . . . . . . . . . . . . . . . . .

32 Atmospheric model

4

2.1

General features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

Discretization and resolution . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.3

Physical parameterizations . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.4

Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.5

Tracer advection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.6

Aerosol indirect eﬀects in large-scale condensation . . . . . . . . . . . . . . .

11

2.7

Cloud and aerosol treatments in radiation scheme . . . . . . . . . . . . . . .

123 Land-surface model

15

3.1

General features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.2

Model structure and resolution. . . . . . . . . . . . . . . . . . . . . . . . .

15

3.3

Processes treated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.4

Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

4 River routing model

17

4.1

General features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

4.2

River routing data and resolution . . . . . . . . . . . . . . . . . . . . . . . .

17

4.3

Treatment of inland depression . . . . . . . . . . . . . . . . . . . . . . . . .

17

5 Ocean model

18

5.1

General features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

5.2

Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

5.3

Model geometry and resolution . . . . . . . . . . . . . . . . . . . . . . . . .

19

5.4

Numerical algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

5.5

Physical parameterization . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

6 Sea ice model

25

6.1

General features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

6.2

Model resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

6.3

Thickness category representation . . . . . . . . . . . . . . . . . . . . . . . .

25

i

6.4

Thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

6.5

Freshwater budget . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

6.6

Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

6.7

Lead parameterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

6.8

Air-sea ﬂuxes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

ii

1

Overview: Structure of MIROC

The Model for Interdisciplinary Research on Climate (MIROC), which is the coupled generalcirculation model used in the K-1 project, consists of ﬁve componet models: atmosphere,
land, river, sea ice, and ocean. The atmospheric component interacts with the land and sea
ice components. The air-sea exchange is realized exclusively between the atmosphere andsea ice components, not directly between the atmosphere and ocean components, and theocean component interacts only with the sea ice component. That is, air-sea ﬂux at ice-free
grids is consequently passed to the ocean component without modiﬁcation, but it is ﬁrst
passed to the sea ice component. The river component receives ground runoﬀ water fromthe land component and drains riverine runoﬀ water into the sea ice component. Lakes aredealt with by the sea ice and ocean components.
Two MIROC setups of diﬀerent resolution are used in the K-1 project and are describedherein. The higher resolution setup is referred to as “HI” and the lower resolution one as“MID” hereafter.

1.1

Flux coupler

Model resolution and grid system are independently chosen for each component, and time
integration also proceeds independently. Note that the sea ice and ocean components sharethe same horizontal resolution and grid system, and their time integration is synchronized.
The ﬂuxes of heat, water, and momentum and the physical quantities necessary for theﬂux calculation (surface air temperature for the air-sea heat ﬂux calculation in the sea ice
component, for example) are passed from one component to another with a prescribed timeinterval. The exchaged ﬂuxes are stocked in the “coupler” and time-averaged over thatinterval, while the latest values are passed at the timing of exchange for the other referenced
variables. The sea ice and ocean components use the same time step, and the exchange isrealized at every time step.
The coupler grid system is generated by dividing each atmospheric grid by meridians andlatitude circles, as the resolution of the land and sea ice components is no lower than thatof the atmospheric component in usual applications. The atmospheric variables necessary
for the surface ﬂux calculation are interpolated onto the coupler grid system. As for theﬂuxes calculated in the atmospheric component and passed to the land or sea ice component(precipitation, for example), the same ﬂux value is assigned to each of the corresponding
coupler grids. The atmosphere-sea ice coupler grid system is also used for the river-sea
ice exchange, and the atmosphere-land coupler grid system is also used for the land-river1

exchange.
Atmosphere-land coupling The atmosphere-land coupler grid system always coincides
with the land component grid system, so the land surface ﬂuxes calculated in the landcompoent model at each grid are simply passed to the atmosphere-land coupler grids, andeach atmospheric grid takes up the ﬂuxes averaged over the corresponding coupler grids.
The number of division for the coupler grid generation is two in both zonal and meridionalfor HI, and is one for MID.
Atmosphere-sea ice coupling The atmosphere-sea ice coupler grid system does notnecessarily coincides with the sea ice component grid system. Each coupler grid generally
covers more than one sea ice component grid, some of which partially and others entirely.
For each coupler grid, therefore, covered sea ice grids and their fractional areal coverage aretabulated beforehand. Likewise, covered coupler grids and their fractional areal coverage
are tabulated for each sea ice model grid. The ﬂuxes and variables are converted between
the coupler grid system and the sea ice model grid system by using these tables and areal
weights. Fluxes are exchaged without any leak through this procedure. The number ofdivision for the coupler grid generation is four in zonal and six in meridional for HI, and is
two both in zonal and meridional for MID.
Land-river-sea ice coupling The river component receives ground runoﬀ water from theatmosphere-land coupler and passes riverine runoﬀ water to the atmosphere-sea ice coupler.
The river component model grid system is generally diﬀerent from either of the two coupler
grid systems, so ground runoﬀ water from more than one atmosphere-land coupler grid ispassed to a single river grid, and riverine runoﬀ water from more than one river grid ispassed to a single atmosphere-sea ice coupler grid. The conversion between the grid systems
is achived by a method similar to the case of the air-sea exchange, and the total water isconserved through this procedure.
Coastal boundary The coastal boundary is deﬁned by the sea ice component grid systemand can be chosen independently of the atmospheric component’s resolution and grid system.
Therefore, land and sea can co-exist in a single atmospheric model grid. Area of coastal land
grids is adjusted so that the sum of the land area and the sea area covered by a single
atmospheric model grid becomes the same as the area of that atmospheric grid.

2

Coupling interval The coupling interval of HI is 3 hours for all the exchanges. The
inteveral for MID is 3 hours for the exchange between sea ice and others, and is 20 minutes
for the other exchanges.

1.2

Coupling strategy for parallel computing

The multi-program multi-data (MPMD) method is used for parallel computing, where thecomponent models are divided into two parts, the atmosphere-land-river part and the sea
ice-ocean part, and these two parts run in parallel on diﬀerent groups of processors. It
means that the coupled model is divided at the atmosphere-sea ice coupler. Necessary data
are exchanged between the two parts at the timing of the atmosphere-sea ice coupling. The
number of processors for each part is so chosen that the computational time required for the
integration over a given period becomes the same between the two parts.

3

2

Atmospheric model

The atmospheric general circulation model used herein is the CCSR/NIES/FRCGC AGCM
version 5.7, which is based on the CCSR/NIES AGCM described in Numaguti et al. (1997;
hereafter, N97). Here, only the outline of the model along with the modiﬁcations done after
N97 was written are described.

2.1

General features

The dynamical part of the CCSR/NIES/FRCGC AGCM solves the primitive equations on
a sphere using a spectral transform method. The vertical coordinate used in this particularversion is σ (pressure normalized by surface pressure), while σ-pressure hybrid coordinate
has also been implemented in the model. For the advection of tracers, a ﬁnite diﬀerence
scheme with the piecewise parabolic method (PPM) is used in conjunction with the ﬂux
form semi-Lagrangian scheme (Lin and Rood, 1996), which is a new feature not having
appeared in N97 and will be described in Subsection 2.5.
The prognostic variables are zonal and meridional velocity, temperature, surface pressure
and mixing ratios of arbitrary number of tracers including water vapor (speciﬁc humidity)
and cloud liquid water. When the interactive aerosol module (sSPRINTARS, a simpliﬁed
version of SPRINTARS: Spectral Radiation-Transport Model for Aerosol Species; Takemura
et al., 2000, 2002) is activated, several aerosol species are also included in the tracers. In
addition, cloud base mass ﬂux of cumulus convection is also prognosed in this particularversion.
The physical parameterizations included are a radiation scheme with two-stream kdistribution method, a prognostic Arakawa-Schubert cumulus scheme, a prognostic cloudwater scheme for large-scale condensation, a turbulence closure scheme with cloud eﬀect,
and an orographic gravity wave drag scheme. What is remarkably new in the physical parameterizations is the treatment of direct and indirect radiative eﬀects of aerosols. The direct
eﬀect and the ﬁrst indirect (cloud albedo) eﬀect is treated in the radiation scheme, while thesecond indirect (cloud lifetime) eﬀect in the large-scale condensation scheme. Those will be
described in Subsections 2.6 and 2.7.

2.2

Discretization and resolution

The spatial discretization method used are a spectral transformation method with Gaussian
grid in horizontal and a grid diﬀerentiation (Arakawa and Suarez, 1983) in vertical.

4For the temporal discretization, a leap-frog scheme is basically used. The integration are
divided into two steps: the dynamics step treating adiabatic advective processes in resolvablescale and the physics step treating the other processes. The time step is variable dependingon maximum wind speed to meet the CFL condition.
In the higher-resolution setup (HI), the horizontal resolution is T106 spectral truncation,
which is approximately equivalent to 1.125◦ longitude-latitude grid. The model has 56
vertical σ layers with relatively ﬁner vertical resolution in the planetary boundary layer
and around the tropopause. The σ values of the vertical layer boundarys (half-levels) are:1.000, 0.995, 0.988, 0.978, 0.966, 0.951, 0.933, 0.912, 0.889, 0.863, 0.835, 0.805, 0.769, 0.725,
0.673, 0.6130, 0.5510, 0.4880, 0.4330, 0.3920, 0.3580, 0.3290, 0.3030, 0.2800, 0.2590, 0.2400,
0.2220, 0.2050, 0.1890, 0.1750, 0.1620, 0.1500, 0.1380, 0.1280, 0.1180, 0.1090, 0.1010, 0.0933,
0.0863, 0.0798, 0.0737, 0.0682, 0.0630, 0.0578, 0.0523, 0.0467, 0.0410, 0.0356, 0.0300, 0.0242,
0.0187, 0.0139, 0.0097, 0.0063, 0.0038, 0.0022 and 0.0000. The height of the model top isapproximately 40km.
In the lower-resolution setup (MID), the horizontal resolution is T42, approximately
equivalent to 2.8125◦. The model has 20 vertical σ layers with relatively ﬁner verticalresolution in the planetary boundary layer. The σ values of the vertical layer boundarys(half-levels) are: 1.000, 0.990, 0.970, 0.930, 0.870, 0.790, 0.700, 0.600, 0.500, 0.410, 0.330,
0.260, 0.200, 0.150, 0.100, 0.070, 0.050, 0.040, 0.030, 0.020 and 0.0001. The height of themodel top is approximately 30km.
Related to the diﬀerence in horizontal and vertical resolution and the model top height,
diﬀerent strength of horizontal diﬀusion is speciﬁed in HI and MID. The e-folding time scale
for the horizontal diﬀusion is 6 hours in HI and 24 hours in MID, while stronger diﬀusion isapplied in the top 4 layers in HI.

2.3

Physical parameterizations

Each of the physical parameterizations adopted are described below in the order in which
they are applied in the model. The prognostic variables are updated three times in eachphysics step due to the adjustment-like nature of moist processes: the ﬁrst is after cumulus
convection, the second after large-scale condensation, and the third is at the end of thephysics step.
Cumulus convection The cumulus parameterization scheme is based on Arakawa andSchubert (1974) with several simpliﬁcations described in N97. The closure assumption is5changed from the diagnostic closure used in N97 to a prognostic closure based on Pan andRandall (1998), in which cloud base mass ﬂux is treated as a prognostic variable. An empirical cumulus suppression condition introduced in Emori et al. (2001) is adopted. Cumulus
convection is suppressed when cloud-mean ambient relative humidity is less than a certain
critical value, which is set to 0.8 in the present version. The momentum transport by cumulus convection based on Gregory et al. (1997) has been implemented and is activatedonly in MID. Since this scheme seems to hinder the organization of tropical precipitation in
HI, it is temporary suppressed in HI. Other modiﬁcations include: an improved downdraft
treatment in which saturated downdraft is more properly represented, gradual melting of
precipitation to avoid too strong temperature kink at melting levels, and a column-based aposteriori water mass ﬁxer to remove negative humidity.
Large-scale condensation The large-scale condensation scheme describes grid-scale condensation and precipitation processes and gives condensational heating, precipitation, cloud
fraction and change of water vapor and cloud liquid water. The scheme has been developed
based on the scheme of Le Treut and Li (1991), in which subgrid probability distribution of
total water mixing ratio in each grid box is assumed as a uniform distribution. The rangeof the uniform distribution was prescribed as a function of height in the version describedin N97 (not explicitly written in N97). In the present version, the range is enhanced wheredetrainment of cumulus convection occurs for an eﬃcient formation of anvil clouds. The
phase of cloud (solid or liquid) is diagnosed according to temperature, T . That is, the ratioof liquid cloud water to total cloud water, fliq , is given by:fliq = min(max((T − Tw )/(Ts − Tw ), 0), 1),
where Ts and Tw are set to 273.15 K and 258.15 K, respectively, in the present version.
However, this diagnosis should be applied only to locally condensed cloud. The fallen ice
cloud should be continued to be treated as ice, and treated as rainfall when it melts. To try
to distinguish the locally condensed cloud and the fallen ice, it is assumed that the ratio of
the former to the total is proportional to the generation tendency of the former to that of
the total. The second indirect eﬀect of aerosols (cloud lifetime eﬀect) is also represented,
which will be described in a separated subsection (Subsection 2.6). Two uncertain parameters primarily used for callibration of the top-of-the-atmosphere radiation budget of theatmospheric model (when realistic SST is prescribed) reside in this part. They are ’vice0’
and ’b1’, and respectively related to ice fall velocity and auto-conversion time scale of liquid
precipitation. In HI and MID, diﬀerent values are speciﬁed for these two parameters.
6

Radiative transfer The radiative transfer scheme employed is based on the two-stream
discrete ordinate method and the k-distribution method. The basic description of the originalradiation code of CCSR/NIES/FRCGC AGCM (hereafter MSTRN) is written in Nakajima
et al. (2000). In the present version, the total number of bands (channels) is 18(37), 11(30)for IR and 7(7) for visible-UV. The number of the band tuned by k-distribution is 13.
The treatments of clouds and aerosols have been improved since Nakajima et al. (2000)
was written especially related to the direct and indirect eﬀects of aerosols. Those will be
described in a separated subsection (Subsection 2.7).
Surface ﬂux The calculation of surface ﬂuxes, including turbulent ﬂuxes of momentum,
sensible and latent heat, and short-wave and long-wave radiative ﬂuxes, has been separated
from the atmospheric model. The atmospheric model receives those ﬂuxes from the ocean
model and the land-surface model through a ﬂux-coupler and uses them as lower boundary
conditions of vertical diﬀusion and radiative transfer.
Vertical diﬀusion The level 2 scheme of turbulence closure model by Mellor and Yamada(1974, 1982) is used for the subgrid vertical ﬂuxes of prognostic variables. The eﬀect offree convection described in N97 has been removed. The treatment of cloud eﬀect for the
turbulent ﬂux is changed from that described in N97 to one based on Smith (1990). In the
new scheme, the thermodynamic variables to be transported by turbulence are total waterand liquid water potential temperature, instead of water vapor and potential temperature
as in the old scheme. After the vertical diﬀusion is calculated, cloud water, water vapor
and temperature are diagnosed in a manner consistent with the large-scale condensation
scheme. A non-local turbulence closure scheme based on Holtslag and Boville (1993) has
also been implemented and can be used in conjunction with the Mellor-Yamada level 2, but
not activated in the present version.
Internal gravity wave drag The eﬀects of orographically exited subgrid scale internal
gravity wave is parameterized following McFarlane (1987). No modiﬁcation has been madeover N97. The strength of this drag is obviously resolution dependent. Diﬀerent parametervalues are used in HI and MID.

2.4

Boundary conditions

The lower boundary conditions required are surface height and the standard deviation of thesurface height over each grid box, which are generated from the USGS GTOPO30 dataset.
7

The three dimensional (or zonal-mean) distribution of ozone concentration is also required.
When the interactive aerosol module is not activated, the three dimensional distributionsof mass concentration of various aerosol species and number concentration of total aerosol
available for cloud condensation nuclei should be prescribed, while various emission data isrequired when the aerosol module is activated.

2.5

Tracer advection

The conservation law in ﬂux form

The ﬂux form of the conservation law for aconservative scalar q in multi–dimensional coordinates is expressed as eq.(2.2) in Lin andRood (1996; hereafter LR96). Discretizing by x = xi (i = 1, 2, 3..), y = yj (j = 1, 2, 3..), andσ = σk (k = 1, 2, 3..), eq.(2.2) in LR96 can be rewritten as∂P S qi,j,k
=
∂t

1
x
(F x 1 − Fi+1,j,k )
2
Δxi,j,k i− 2 ,j,k
1y
(F y 1 − Fi,j+
+
1 )
,k
2
Δyi,j,k i,j− 2 ,k
1σ
(F σ 1 − Fi,j,k++
1 ).
2
Δσi,j,k i,j,k− 2(1)xHere, P S denotes the surface pressure, Fi−
denotes the mass ﬂux of tracer constituents for
1
,j,k
2

the x–direction on the boundary between cell (i, j, k) and cell (i − 1, j, k) (cf. Fig.1 in LR96),
and Δxi,j,k denotes the width of cell (i, j, k) for the x–direction. The total constituent mass
(global integral of constituent density) is generally conserved in an operator–split ﬂux–form
scheme without a posteriori restoration. The precision and characteristics (e.g., stability andmonotonicity) of this advection scheme depends on the selection of computational schemex
for the calculation of mass transfer across the cell (Fi−). CCSR/NIES/FRCGC AGCM
1
,j,k
2

x
can choose between two types of van Leer–type schemes for the calculation of Fi−
, by
1
,j,k
2

considering computational cost and the precision. One is the 4th–order monotonic van Leerscheme (called as FFSL–2 in LR96) and the other is the Piecewise Parabolic Method (PPM)
scheme with monotonic constraint (FFSL–3 in LR96). Lin et al. (1994) describes generalformulation of the 4th–order monotonic van Leer scheme, and Carpenter et al. (1990)
and Colella and Woodward (1984) describe general formulation of the PPM scheme. The
monotonicity constraint for the PPM scheme is described in Appendix C of LR96. Each ofthese schemes
1. compute ﬂuxes based on the sub-grid distribution in the upwind direction
2. generate no new maxima or minima (i.e., monotonicity is satisﬁed)
8

3. generate no negative value
4. are stable when the Courant number in each direction is less than 1.
The PPM scheme needs more computational cost compared to the 4th–order van Leerx
scheme, but it can calculate Fi−
with less computational dissipation. Therefore, the
1
,j,k
2

PPM scheme is applied for the k–1 coupled–model calculation.
Application to a spectral model For the calculation of the time–averaged mass ﬂux
across the cell boundary, the winds and the tracer distributions are staggered in the Arakawa
C–grid (cf., Mesinger and Arakawa, 1976). The horizontal winds at the cell boundary, ui− 1 ,j,k2

and vi,j− 1 ,k , are reconstructed by using the mass convergence ﬁeld in the spectral model and2

the discretized continuity equation:
S∂Pi,j,k
1
ϕ
ϕ
λ
λ
σ
σ
=
{(Vi−
− Vi+
) + (Vi,j−− Vi,j+1
1
11 ) + (Vi,j,k− 1 − Vi,j,k+ 1 )}.
,j,k
,j,k
,k
,k
2
2
2
2
2
2
∂t
ΔDj,k

(2)

ϕ
λ
σ
, Vi,j−denote zonal, meridional, and vertical mass–weighted
Here, Vi−1
1 , and V
i,j,k− 1
,j,k
,k
2

2

2

wind at the cell boundary, respectively. That is,
λ
S
1
Vi−
1
,j,k = (P u)i− ,j,k Δyj Δσk(3)ϕ= (P S v)i,j− 1 ,k Δxj− 1 Δσk
Vi,j−
1
,k

(4)

σ
S
Vi,j,k−1 = (P σ̇)i,j,k− 1 Δxj Δyj .

(5)

2

2

2

2

2

2

2

ΔDj,k denotes the cell volume, and Δxj , Δyj , and Δσk denote zonal, meridional, vertical
width of the cell, respectively. That is, ΔDj,k = a cos ϕj ΔλΔϕj Δσ, Δxj = a cos ϕj Δλ, andΔyj = aΔϕj .
The following are the procedure for the calculation of tracer advection in the staggering–
grided horizontal and vertical wind ﬁelds:
1. Surface pressure P S and horizontal wind v at time step t + Δt are predicted in the
spectral model.
2. The horizontal component of mass ﬂux divergence at time step t is calculated by using
spherical harmonics. The mass ﬂuxes at time step t are reconstructed from the values
at t + Δt and t − Δt, because CCSR/NIES AGCM applies semi–implicit scheme forthe time–integration of surface pressure. Zonal and meridional component of mass ﬂux
divergence are:Cx = −∂∂
1
1(P S u), C y = −(P S v cos ϕ).
a cos ϕ ∂λa cos ϕ ∂ϕ9

3. By applying Cx and Cy into eq.(3) and eq.(4), V λ , V ϕ are calculated by using:
λ
λ
x
Vi−
1
,j,k − Vi+ 1 ,j,k = Ci,j,k ΔDj,k

(6)

ϕ
ϕ
y
Vi,j−
− Vi,j+= Ci,j,k
ΔDj,k .
1
1
,k
,k

(7)

2

2

2

2

The boundary conditions for eq.(6) and eq.(7) are V ϕ = 0 at the North Pole and theSouth Pole, and zonal average of V λ is equal to that in the spectral model grid.
By integrating eq.(2) from the top to the bottom of the vertical σ domain,S∂Pi,j,kxx
ΔDj,k = (Ci,j,k+ Ci,j,k
).
∂t
k
k

(8)

V σ is calculated from eq.(2) and eq.(8), with the boundary condition of V σ = 0 at thesurface and the top of the atmosphere.
4. The mass ﬂux for tracer constituents F x , F y , and F σ are calculated by using van
Leer–type scheme (4th order monotonic van Leer scheme or monotonic PPM scheme)from the mass–weighted winds V λ , V ϕ , and V σ . If the absolute value of the Courant
number exceeds 1, a Flux–form semi–Lagrangian (FFSL) treatment is applied, i.e., themass ﬂux is separately calculated as the integer ﬂux and the fractional ﬂux, following
to LR96. Since the meridional and vertical ﬂux are negligibly small compared to thezonal ﬂux, FFSL treatment is considered only for F x at present.
5. P S q at time step t + Δt is calculated from eq.(1). qt+Δt can be obtained by dividing(P S q)t+Δt by the surface pressure at time step t + Δt. There is small quantity of diﬀerS
S
ence between Pt+Δt
from eq.(8) and Pt+Δt
in the spectral model, because CCSR/NIESS
AGCM applies the semi–implicit time integration scheme. Pt+Δt
from eq.(8) is appliedat present for the consistency of mass advection. Mass conservation is not strictly
satisﬁed because of the discrepancy between the surface pressure in the spectral modelS
and Pt+Δtfrom eq.(8).

The mass ﬂuxes into/out of polar caps

The mass ﬂuxes into/out of polar caps are cal-culated by using the semi–Lagrangian scheme in the polar stereo projection (cf. Figure 2.5).
The horizontal average at the highest latitude band is assumed to be preserved before/after
ﬂux calculation for the mass conservation. The sequence of calculation is:
1. Zonal average of P S q at time step t is calculated at the highest latitude band, and isassumed to equal P S q at the pole.
10y:λ=90u(i,j)

Grid Pt.(i,j=1)

Upstream Pt.
(xs,ys)v(i,j)

N.P

x:λ=0
YP

YP*R2

Figure 1: Conceptual ﬁgure for the ﬂux on pole-most grids.

2. Horizontal wind at the highest latitude bands is projected into the orthogonal coordinate system centering around the pole, and q at time step t + Δt is estimated by using
the value at the “departure point”.
3. Zonal average of (P S q)t+Δt is ﬁxed to that of (P S q)t .

2.6

Aerosol indirect eﬀects in large-scale condensation

The present version considers both the ﬁrst and second aerosol indirect eﬀects only forstratus water cloud. The cloud droplet number concentration Nc is diagnosed referring tothe observation by Martin et al. (1994):
Nc =Na Nm,Na + Nm

where Na is the aerosol particle number concentration excluding soil dust and black carbon,
and  and Nm are constants. The monthly mean data of the three-dimensional aerosol
number concentration used here is prescribed oﬀ-line by SPRINTARS, or they are calculated
on-line at each time step when the interactive aerosol module is activated. If the function of
the cloud droplet size distribution is ﬁxed, its eﬀective radius ref f is decided as follows:3 ρlref f = k
4πρw Nc
11

1
3

,

where ρ is the air density, l is the cloud water mixing ratio, ρw is the water density, and k
is an empirical constant. The Berry’s parameterization is adapted to the precipitation rate
P , which depends on not only l but also Nc (Berry, 1967):
P =−

dl
αρl2
,
=
dt
β + γ Nρlc

where t is time and α, β, and γ are constants.

2.7

Cloud and aerosol treatments in radiation scheme

The optical parameters for seven species of particulate matter, two types of cloud (water
cloud and ice cloud) and ﬁve aerosol species (soil dust, black carbon (BC), organic carbon(OC), sulfate, and sea salt) are included. The classiﬁcation of aerosols is based on SPRINTARS (Takemura et al., 2000, 2002). The parameters are calculated oﬀ-line, and the GCM
access them through a table.
Aerosol microphysics and its radiative eﬀect on water cloud As mentioned in the
Subsection 2.6, the radius of water cloud varies depending on the aerosol particle number
concentration. The monthly mean data of the three-dimensional aerosol mass and number
concentrations used in the radiative process are prescribed oﬀ-line by SPRINTARS, or they
are calculated on-line at each time step when the interactive aerosol module is activated.
The hygroscopic growth of aerosols is expressed according to Tang and Munkelwitz (1994)
for sulfate and Hobbs et al. (1997) for OC, and they are assumed to be homogeneouslyintermingled with water. The relationship between the relative humidity and its volume
mode radius is shown in Table 1.
Table 1: volume mode radius rg of OC and sulfate depending on relative humidity (RH).
RH (%)0rg (μm) (OC)0.237
rg (μm) (sulfate) 0.243

50
0.256
0.298

70
0.261
0.333

80
90
0.342 0.401
0.361 0.427

95
0.456
0.550

98
0.650
0.683

99
0.740
0.809

Optical parameters for clouds and aerosols The refractive indices of water and dry
aerosols used in preparing the optical parameter table are according to d’Almeida et al.
(1991) and WCP-55 (1983), but their imaginary part of soil dust at visible wavelength isone-fourth of WCP-55 (1983) in this model, which is ﬁtted to recent observations. The
12

coeﬃcients of extinction and absorption and the ﬁrst ﬁve orders of Legendre moments ofthe volume scattering phase function in each wavelength band are written in the table. We
assume that the shape of the particulate matters spherical. The distribution of particle radiiis calculated assuming a lognormal distribution and taking the values shown in Table 2 as
the volume mode radius and standard deviation of the distribution. The optical parameters
are calculated for the eight sizes of volume mode radius for each species of particulatematter. For the hygroscopic aerosols (OC and sulfate aerosols), the optical parameters
are calculated for eight radii in accordance with the relative humidity taking account of
hygroscopic particle growth as described in the previous section (Table 1), and then they
are interpolated depending on the relative humidity in each model grid. The eﬀect of theparticle size on the radiative ﬂux is expressed by the interpolation of the moments of thevolume scattering phase function of the referenced volume mode radius in the table.
Table 2: volume mode radius rg and standard deviation σ for each particulate matter. ∗ see
Table 1.
rg (μm)σwater cloud ice cloudvariable
40.0
1.5
1.5

soil dust
BC
4.00
0.03902.5
2.0

OCvariable∗
1.8

sulfatevariable∗
2.0

sea salt
2.00
2.5

Vertical cloud overlapping In the original MSTRN, the semi-random method was used
to treat cloud overlap (Nakajima et al.,2000). While this computationally eﬃcient method
works well for an AGCM, the appropriateness of random overlap assumption depends onvertical resolution of the model.
This is because, in the random overlap assumption, the ratio of cloud overlapping increases as the number of vertical layers increases. When a mass of cloud exist in the continuous layers, it is not reasonable. When the model resolves a cloud, the adjacent layers
which contain cloud should be treated assuming maximum overlap.
Several diﬀerent algorithms have been proposed which compute the maximum-random
overlap assumption, in which adjacent cloud layers overlap maximally and discrete cloudlayers are independent. We adopted the method of Geleyn and Hollingsworth (1979, hereafterGH). This method allots the radiative ﬂuxes to clear and cloudy segments according to theratio of the cloud fractions in adjacent layers. As a result, the radiative ﬂuxes are obtained
under the condition of maximum-random overlap. Since we calculate the multiple scatteringof light between model layers by the adding method (Nakajima et al., 2000), we cannot apply
13the original GH method to MSTRN directly; therefore we re-formulated the GH method.
The new formulation is mathematically equivalent to the original, and is similar to that of
the four-stream approximation (Shibata and Uchiyama, 1992).

14

3

Land-surface model

The land-surface model used herein is MATSIRO (Minimal Advanced Treatments of Surface
Interaction and RunOﬀ), which is described in Takata et al. (2003).

3.1

General features

MATSIRO represents energy and water exchange between land-surface and atmosphere.
MATSIRO receives temperature, speciﬁc humidity, wind speed, and pressure of the lowest atmospheric layer, and precipitation and downward short-wave and long-wave radiation
ﬂuxes from the atmospheric model, while it provides turbulent ﬂuxes of momentum, sensible
and latent heat, and upward short-wave and long-wave radiation ﬂuxes to the atmosphericmodel. It also provides runoﬀ ﬂux to the river routing model. The prognostic variables of
MATSIRO are canopy temperature and surface temperature of snow-free and snow-covered
areas, canopy water content, snow amount, snow temperature, snow albedo, soil temperature, soil moisture content, and frozen soil moisture content.

3.2

Model structure and resolution

The model consists of two parts: ﬂux part and integration part. The ﬂux part is invokedby the atmospheric model in each time step of the atmospheric model. It provides ﬂuxes tothe atmosphere and updates canopy and surface temperature. The integration part updates
canopy water content and snow and soil variables. Its time step is 3 hours.
The model has 1 canopy layer, 5 soil layers, and variable number of snow layers (0–3
depending on the snow amount). The thickness of the soil layers from the top are: 5cm,
20cm, 75cm, 100cm and 200cm. Canopy and surface skin layer do not have heat capacity.
In HI, the horizontal grid structure of the model is 2×2 tiling of each atmospheric model
grid box, that is, the resolution is approximately 0.5625◦ in longitude and latitude. In MID,the horizontal resolution and the grid structure are the same as those of the atmosphericmodel, that is, T42 (∼2.8125◦ ).

3.3

Processes treatedSurface ﬂux The turbulent ﬂuxes from canopy and surface are calculated based on Watanabe (1994). The ﬂuxes over snow-free and snow-covered areas are calculated separately andaveraged afterward weighted with area fractions. The stomatal conductance is evaluated on
the basis of a photosynthesis-stomatal conductance model (e.g., Farquhar et al., 1980; Ball
15

1998).
Canopy water The canopy water content is determined by the balance of canopy interception of precipitation, canopy evaporation and dripping.
Snow The snow amount is determined by the balance of snow fall, snow melt and sublimation. The snow temperature is calculated by a heat conduction equation. The snow albedo
calculation is based on Wiscombe and Warren (1980).
Runoﬀ The TOP model approach (Beven and Kirkby, 1979) is adopted with simpliﬁcation.
Saturation excess runoﬀ, inﬁltration excess runoﬀ, and base ﬂow are evaluated separately.
Soil Heat conduction, water ﬂow due to gravity and suction, freezing and thawing of soil
water and root uptake of soil water are considered.

3.4

Boundary conditions

The boundary conditions required are land-cover type, soil type, leaf area index (LAI),
soil albedo for visible and near-IR wavelength bands, mean surface slope, and the standard
deviation of surface height. Most of them are generated from the ISLSCP Initiative I dataset,
while the land-cover type is based on USGS GLCC and the standard deviation of surface
height is based on USGS GTOPO30. The LAI should be prescribed by a monthly dataset,
while the others are temporary constant.

16

4River routing model

4.1

General features

The river routing model transports runoﬀ water from the land-surface model to the ocean
model at river mouths. The scheme is based on Miller et al. (1994), which prognose riverwater storage in each grid box. A globally uniform and temporary constant river ﬂow speedof 0.3 m/s is used for simplicity. The ﬂow of glacier (the snow amount exceeding a certain
critical value is treated as glacier in the land-surface model) is also treated in this model,
though the freezing and thawing of river water is not considered.

4.2

River routing data and resolution

In HI, the river routing dataset used is TRIP (Total Runoﬀ Integrating Pathways; Oki andSud, 1998). It speciﬁes river ﬂow direction as one of the 8 directions (N, NE, E, SE, S,
SW, W, NW) at each grid box. The resolution of the model and the routing dataset is0.5◦ in longitude and latitude. In MID, the resolution of the model is the same as that of
atmospheric model, that is, T42 (∼2.8125◦ ). The river ﬂow direction dataset, similar in
structure to TRIP, is the one developed by Dr. Shinjiro Kanae (personal communication).
The time step is 3 hours for both HI and MID.

4.3

Treatment of inland depression

In the TRIP dataset, there are some inland points which have no outlets (inland depression).
In HI, which uses TRIP, since evaporation of river water or inﬁltration of river water into
soil is not considered, river water tends to be accumulated at those points. As a temporary
expedient to avoid this accumulation, the river water in the inland depression points issprinkled over surrounding grids. There is no inland depression in MID.

17

5

Ocean model

The ocean general circulation model used herein is CCSR Ocean Component Model (COCO)
version 3.4. More detailed description of COCO can be found in Hasumi (2000).

5.1

General features

The basic equations of COCO is the primitive equations on a sphere, where the Boussinesq
and hydrostatic approximations are adopted. The model explicitly represents the sea surface
elevation. The vertical coordinate of the model is a hybrid of σ (normalized depth) and z(geopotential height): the former is applied between the sea surface and a ﬁxed level inthe upper ocean, and the latter below. The σ-coordinate is introduced in order to avoid
disappearance of the top level. The model is formulated horizontally on the spherical coordinate system, and the rotation of the model grid system is available to avoid the coordinate
convergence in the Arctic Ocean.
In solving the momentum equations, the horizontal velocity is decomposed into thebarotropic mode (vertically integrated velocity) and the baroclinic mode (deviation fromthe vertical mean) in order to take a separate care for high-speed external gravity waves. A
shorter time step is applied for the barotropic mode equations than for the remaining ones
in numerical time-integration.

5.2

Discretization

The model equations are ﬁnite-diﬀerenced in both time and space. The grid spacing and thetime step are ﬁxed in time.
The model’s prognostic variables are arranged horizontally on the Arakawa B staggard
grid system. Coast lines are deﬁned by the lines connecting velocity points. The partial step
formulation (Adcroft et al., 1997) is applied for bottom topography representation, whichallows to set the depth of the sea ﬂoor irrespective of the model vertical grid system.
Numerical time integration of the model’s prognostic equations other than for the barotropic
mode are basically achieved by the leap-frog and Euler-backward schemes: nine steps of thelatter is applied after every ninety steps of the former, and one step of the Euler-forward
scheme is used after the application of the Euler-backward scheme to prepare the valuesof the variables at two time levels. The tracer advection term is, however, integrated byuse of the Euler-forward scheme, as an upstream-weighted algorithm is adopted. The horizontal diﬀusion and viscosity terms are also integrated by the Euler-forward scheme, while18the vertical diﬀusion and viscosity terms are integrated implicitly in time. A semi-implicit
treatment is applied for the Coriolis term. For the barotropic mode equations, only the
Euler-backward scheme is used as proposed by Killworth et al. (1991).

5.3

Model geometry and resolution

HI-resolution The model spherical coordinate system is rotated so that the singularpoints of the coordinate are located on land. The north pole of the model coordinate system
is shifted at 50◦ W, 77◦ N on Greenland. The horizontal resolution is 0.28125◦ zonally and0.1875◦ meridionally, both in the rotated coordinate system. There are 47 vertical levels,
excluding the bottom boundary layer, and 8 of which are within the σ-coordinate. The
vertical grid spacing varies with depth, whose values are 2.5, 2.5, 3, 4, 5, 6, 7, 8, 10, 12, 14,
16, 18, 20, 25, 30, 35, 40, 50, 60, 80, 100, 120, 140, 160, 180, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 200, 200, 200, 250, 250, 250, 250, 250, and 250 m, from the
top to the bottom. Note that the values for the σ-coordinate region are under the zero sea
level elevation. The time step of the model integration is 3 seconds for the barotropic mode
equations and 3 minutes for the others.
The model bathymetry is generated from a widely used 2-minute resolution bathymetry
data (ETOPO2). The depths of the ocean ﬂoor are ﬁrst interpolated on the model grids,
and then are spatially smoothed at high latitudes. The model geometry in the rotated grid
system is shown in Figure 2. Most of the narrow straits connecting major basins and inlandseas are naturally resolved at the current resolution. However, the Straits of Gibraltar is
artiﬁcially widened. The depth of the Strait of Gibraltar is set at 668 m (23rd level).
MID-resolution The model coordinate system coincides with the geographical sphericalcoordinate system. The zonal resolution is 1.40625◦ (= 360/256). The meridional resolution
is about 0.56◦ at latitudes lower than 8◦ , 1.4◦ at latitudes higher than 65◦ , and smoothly
changes in between. There are 43 vertical levels, excluding the bottom boundary layer, and
8 of which are within the σ-coordinate. The vertical grid spacing varies with depth, whose
values are 5, 5, 5, 5, 5, 5, 7.5, 7.5, 10, 10, 10, 10, 15, 20, 25, 25, 30, 50, 75, 75, 100, 100, 150,
150, 150, 150, 200, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250,
250 m, from the top to the bottom. Note that the values for the σ-coordinate region are
under the zero sea level elevation. The time step of the model integration is 30 seconds forthe barotropic mode equations and 20 minutes for the others. Sea surface elevetion, tracers,
and velocity components are ﬁltered in the zonal direction in the Arctic Ocean to use a time

19

Figure 2: Model ocean bathymetry for HI in the rotated coordinate system.
step longer than the limit of the CFL condition in the Arctic Ocean. The ﬁltering is done byrepeated operation of the simplest symmetric ﬁnit impulse response ﬁlter, where the ﬁlteredvalue at a certain grid is obtained by the weighted average of the value at that grid and two
adjacent grids with the weight of 2 for the central grid and 1 for the adjacent grids. Care istaken to conserve water volume and tracers.
The model bathymetry is generated from a widely used 5-minute resolution bathymetry
data (ETOPO5). The depths of the ocean ﬂoor are ﬁrst interpolated on the model grids,
and then are spatially smoothed. The model geometry in the rotated grid system is shownin Figure 3. The Bering Strait is represented by a two-grid gap, which means there is only
one velocity grid point at the strait. The water pathway through the Canadian Archipelago
is also represented by artiﬁcially excavating a channel. The Mediterranean, Black, Caspian,and Red Seas and the Hudson Bay are represented as isolated basins. At the straits connecting these basins with the adjacent seas and oceans, sea surface elevation and tracers are
artiﬁcially exchanged by two-way linear damping. The time constant of the damping is 100
days. The depth to which the damping is applied is 1650 m for the Strait of Gibraltar and
45 m for the other straits.

20

Figure 3: Model ocean bathymetry for MID.

5.4

Numerical algorithm

Barotropic mode equations The barotropic mode equations are made up of the vertically integrated momentum equations and the vertically integrated continuity equation.
The formulation follows the method of Killworth et al. (1991). When the barotropic velocity components are combined with the baroclinic ones, an average over two baroclinic
time steps is taken in order to satisfy the continuity equation, as indicated by Dukowicz andSmith (1994). In order to eliminate checker-board noise in the sea surface elevation ﬁeld,
weak horizontal diﬀusion is applied to sea surface elevation, with a care to conserve tracer
quantities in the σ-coordinate region. The coeﬃcient for the diﬀusion is 10 m2 /s for both
HI and MID.
Momentum advection The model employs the momentum advection algorithm of Ishizaki
and Motoi (1999) for HI, which is a pseudo-enstrophy preserving scheme with a consideration for up-/down-sloping advection. The momentum advection algorithm for MID is thesecond-order, energy-preserving one.
Tracer advection The model’s tracer advection is based on the uniformly third-orderpolynomial interpolation algorithm (UTOPIA; Leonard et al., 1993). In its application to
the model, however, the horizontal dimensions and the vertical dimension are separately
21

dealt with, so the vertical algorithm reduces to the quadratic upstream interpolation forconvective kinematics with estimated streaming terms (QUICKEST; Leonard, 1979). The
neglected vertical-horizontal cross term is generally small for hydrostatic models where thegrid aspect ratio (vertical to horizontal) is much smaller than unity.
In the horizontal application of UTOPIA, the ﬂux-integral method (Leonard et al. 1994)
is used to guarantee the tracer quantity conservation. A ﬂux limiter is required for numerical stability of QUICKEST and UTOPIA. There are single-dimensional ﬂux limiters (e.g.,
Leonard, 1991) and multi-dimensional ones. The former has an advantage in terms of positivity preservation, while the latter has an advantage in terms of shape preservation, relative
to the other. A single-dimensional one is adopted herein, as the multidimensional ones devised by Leonard et al. (1993) resulted in a problem where there is a strong cross-gradientcurrent and its orientation is diagonal to the horizontal grid. Such a feature is often observedin the bottom boundary layer.

5.5

Physical parameterization

Vertical convection Since the hydrostacy is assumed, the vertical convection process
should be appropriately parameterized to remove unstable vertical stratiﬁcation. The way
adopted here is the simplest convective adjustment, where a statically unstable water column
is homogenized instantaneously.
Vertical mixing The vertical mixing of momentum and tracers is represented by a harmonic form. The diﬀusion and viscosity coeﬃcients are calculated by the parameterizationof Noh and Kim (1999), which is based on the level 2.5 turbulence closure of Mellor and
Yamada (1982). The formulation is slightly modiﬁed to include dependence of the turbulent
Prandtl number on the turbulent Richardson number:
P r = P r0 + P r1 Ri,(9)where P r is the turbulent Prandtl number, P r0 is its value at the neutral stability, Ri is
the turbulent Richardson number, and P r1 is the parameter determining the dependence on
Ri. P r0 is taken to be 0.8 as suggested by many experimental studies. P r1 is taken to be 7,
following Kondo et al. (1978). An upper bound is imposed on P r, which is set to 20. The
value of the parameter α, which determines the model’s behavior at high stability, is taken
to be 3. Other model parameters are the same as given by Noh and Kim (1999).
As background diﬀusivity, a minimum value is prescribed for each level. The backgroundviscosity coeﬃcient is 1.0 × 10−4 m2 /s at all the levels. The background diﬀusion coeﬃcient22

varies with depth, from 0.1 × 10−4 m2 /s for the top to 3.0 × 10−4 m2 /s for the deepest
model level. This vertical proﬁle is suggested by Tsujino et al. (2000) to simulate the Paciﬁc
deep overturning circulation with suﬃcient intensity. Latitudinal dependence is assumed
for the background vertical diﬀusivity. Gregg et al. (2003) suggest that the dissipation
rate at low latitudes is signiﬁcantly reduced compared with that for the same internal wave
characteristics at higher latitudes. When it is referenced to the value derived from the Garretand Munk’s (1975) internal wave spectrum at 30◦ , the factor to be multiplied by is modeledby them as

f cosh−1 (N/f ),f30◦ cosh−1 (N0 /f30◦ )(10)

where f is the Coriolis parameter, f30◦ is its value at 30◦ , N is the buoyancy frequency, andN0 (= 5.24 × 10−3 /s) is the reference stratiﬁcation assumed in the model of Garret andMunk (1975). This functional form is applied directly to the background diﬀusivity, not tothe dissipation rate, in our model. Since N is not known a priori, it is taken to be the sameas N0 for the current application.
Lateral mixing of momentum The Smagorinsky’s (1963) biharmonic viscosity is applied
for the lateral momentum mixing in HI. Its coeﬃcient is dependent on the grid width and thestrain rate, and is controlled by a single nondimensional parameter whose value is taken to
be 4. The harmonic viscosity is used in MID with the coeﬃcient value of 3.0 × 104 cm2 /s. At
latitudes higher than 60◦ , the coeﬃcient is linearly reduced to zero with increasing latutide.
Lateral mixing of tracer The constant coeﬃcient biharmonic diﬀusion is applied as the
background in HI: the value of the coeﬃcient is 1.0 × 109 m4 /s. Harmonic isopycnal diﬀusion
is also applied with the coeﬃcient value of 1.0 × 103 m2 /s.
Since the HI-resolution model’s horizontal grid spacing is not ﬁne enough to resolve
mesoscale eddies at high latitudes, it is necessary to introduce some parameterization foreddy-induced transport of tracers. The harmonic horizontal diﬀusion of isopycnal layer
thickness (Gent et al., 1995) is applied at high latitudes. The value of the thickness diﬀusion
coeﬃcient is 7.0 × 102 m2 /s for latitudes higher than 50◦ , is 0 for latitudes lower than 40◦ ,and linearly changes in between.
For MID, the harmonic horizontal diﬀusion, the harmonic isopycnal diﬀusion, and theharmonic hirozontal diﬀusion of isopycnal layer thickness are applied, whose coeﬃcients are
1.0 × 102 m2 /s, 1.0 × 103 m2 /s, and 7.0 × 102 m2 /s, respectively.

23

Bottom boundary layer Down-sloping ﬂow of dense water, which is a crucial process
for the deep overturning circulation, is not naturally represented at the current model resolution. The bottom boundary layer parameterization of Nakano and Suginohara (2002) isincorporated into the model to represent it. In this parameterization, a bottom boundary
layer of ﬁxed depth is assumed, entrainment at the top of the layer is represented by enhanced vertical diﬀusivity, and the ageostrophic eﬀect of eddies is represented by the bottom
Rayleigh drag. The layer thickness is 80 m for HI and 100 m for MID, the vertical diﬀusionand viscosity coeﬃcients at the top of the layer are 5.0 × 10−4 m2 /s and 1.0 × 10−4 m2 /s,
respectively. The Rayleigh drag coeﬃcient is taken to be the same as the Coriolis parameter
above 2000 m and zero below. Horizontal viscosity is considered within the layer, whose
formulation and coeﬃcient are the same as for the upper layers. No horizontal diﬀusion of
tracers is applied. The bottom boundary layer is applied only at high latitudes, to the northof 49◦ N and to the south of 54◦ S, following Nakano and Suginohara (2002).
Vertical penetration of shortwave radiative ﬂux The shortwave radiative ﬂux absorbed by the ocean penetrates into deeper levels. The attenuation of the shortwave ﬂux inthe ocean is calculated by following the formulation of Rosati and Miyakoda (1988).

24

6

Sea ice model

6.1

General features

Sea ice is treated as a two-dimensional continuum. In each horizontal grid, concentration,mean thickness (over the ice-covered portion), and horizontal velocity components are predicted. Sea ice concentration and thickness change due to thermodynamic growth/decay,
advection, and mechanical deformation of ice ﬂoes.

6.2

Model resolution

The horizontal resolution of the sea ice model is taken to be the same as that of the ocean
model. The time step of the integration is also the same as for the ocean model (for other thanbarotropic mode equations). The time splitting method is used in solving the momentum
equation with the elastic-viscous-plastic rheology (described later), where the time step ischosen to be 3 seconds for HI and 30 seconds for MID.

6.3

Thickness category representation

This sea ice model is classiﬁed as two-category in terms of the representation of subgrid-scalesea ice thickness distribution. In typical two-category models, a minimum value is prescribedfor mean thickness, and this minimum thickness is used as a parameter distinguishing thinand thick ice. Mechanical characteristics are very diﬀerent between thin ice (such as pancakes) and thick ice (such as multi-year ice and ridged ice). In the early stage of new ice
formation, sea ice is often formed as frazil, which does not exhibit such a solid form as an ice
ﬂoe. Even after sea ice takes a solid form such as a pancake, collision of thin ice ﬂoes often
results in rafting, not ridging. From the dynamical perspective, therefore, it is reasonable
to treat ice thinner than a certain threshold categorically as identical to open water. In themodel, newly formed thin ice is instantaneously converted into ice of the minimum or larger
thickness. The minimum thickness in this model is 0.1 m.

6.4

Thermodynamics

The simplest zero-layer thermodynamics (Semtner, 1976) is assumed in the sea ice model,
where the heat content of sea ice is neglected and only the latent heat of melt is taken into
account in the heat budget consideration. Snow cover on sea ice is treated likewise.

25

The vertical heat transfer within sea ice is calculated by
To − Ts(11)
FI = kI ks
kI hs + ks hI
where FI is the vertical heat ﬂux through sea ice (positive upward), hI is the sea ice thickness,
hs is the snow cover thickness, and kI (= 2.04 W/m/K) and ks (= 0.31 W/m/K) are theheat conductivity of sea ice and snow, respectively. To is the ice bottom temperature, whichcoincides with the freezing point of the seawater below, and Ts is the temperature at the
air-ice (snow when exists) interface. The freezing point of seawater, Tf , is given by a functionof salinity as Tf = −0.0543S, where Tf is measured by the degree Celsius and the unit ofsalinity S is the practical salinity unit (psu).
The air-ice/snow interface temperature is determined by the heat balance at the interface.
If the diagnosed air-ice/snow interface temperature goes beyond the freezing point of purewater, the temperature is reset at the freezing point and the heat ﬂux imbalance is consumed
by melting ice/snow.
Sea ice exists only when and where the ocean model’s top level is at the freezing point. If
the seawater temperature under sea ice exceeds the freezing point in the course of the model
integration, the sea ice amount is reduced until the seawater temperature settles down at
the freezing point.
When open water of a partially ice-covered grid is cooled, newly formed thin ice eventually
reduces the mean sea ice thickness. This process can be represented by setting the predicted
sea ice concentration at unity. In a two-category model, however, such a treatment could
conﬂict with the interpretation of thin ice in the model. Still, some increase in concentration
should happen under cooling, even by reducing mean thickness. From this standpoint, it
should be parameterized to what degree concentration should increase under the coolingof open water. When open water of a partially ice-covered grid is warmed, on the other
hand, absorbed water is to be consumed to melt sea ice both from the ice edge (reduces theconcentration) and from the ice bottom (reduced the thickness): the latter due to the heattransfer within the sea surface mixed layer. Thus, in this case, the partitioning between thelateral and basal melting should be parameterized. These are parameterized by following themethod of Mellor and Kantha (1989), with the same parameter values as described therein.

6.5

Freshwater budget

Sea ice is assumed to contain constant salinity, which is taken to be 5 psu in the current
modeling. Melt ponds are not taken into consideration in the model, so the melt water
immediately runs oﬀ into the ocean.
26

The level of the snow-ice interface relative to the sea surface is determined by the densityof snow and ice. When it goes below the sea surface, snow is converted to sea ice until theinterface level coincides with the sea surface (snow-ice formation). In this process, salinity isadded to the newly formed snow-ice to guarantee the salt conservation. It is more physicallyreasonable to reduce the salinity of sea ice, but such a treatment requires to deal with thesea ice salinity as a prognostic variable.

6.6

Dynamics

Sea ice concentration and mean thickness are advected by the sea ice horizontal velocity withconserving the volume. The advection term is numerically integrated by using the weighted
upstream algorithm, which is a mixture of the upstream and centered-in-space algorithms.
Harmonic and biharmonic diﬀusion is taken into account for the sake of numerical stability,
whose coeﬃcients are 1.0 ×102 m2 /s and 1.0 ×108 m4 /s, respectively, for HI. The coeﬃcientsare 1.0 × 104 m2 /s and 1.0 × 1012 m4 /s for MID.
The momentum equation of sea ice is composed of the advection term, the Coriolis
term, the acceleration term due to sloping of the sea surface, the internal stress term, andthe external forcing term (e.g., Hibler, 1979). The external forcing is the sum of windstress and ice-ocean drag. A quadratic form is assumed for the ice-ocean drag term, as in
Hibler (1979), and the drag coeﬃcient is taken to be 0.005. Sea ice concentration is notallowed to exceed unity. When ice volume convergence occurs in a compactly packed grid,
ice thickness increases without changing concentration. This represents the sea ice ridging
process. The mechanical resistance to the ridging is represented by the internal stress term in
the momentum equation, and is formulated by the elastic-viscous-plastic rheology of Hunkeand Dukowicz (1997). The parameters for the internal stress are also taken from Hunke and
Dukowicz (1997), except for the parameter representing the strength of the ice of the unit
thickness and concentration (P ∗ in their paper), which is taken to be 2.0 × 102 N/m.

6.7

Lead parameterization

Leads are generally underestimated in large-scale models due, for example, to the lack ofshort-term forcing variations. Observations indicate that 1–2 % of open water exists evenfor the wintertime central Arctic, and the air-sea heat exchange through such leads accounts
for more than half of the wintertime basin-scale heat budget (e.g., Maykut, 1982).
In this model, lead formation is parameterized by imposing a maximum value for sea
ice concentration. The maximum concentration depends on thickness: it is 0.95 for the27minimum thickness, is unity for ice thicker than 6 m, and linearly changes in between.

6.8

Air-sea ﬂuxes

The sea ice model also functions as the coupled model’s air-sea interface, and the ﬂuxes
exchanged between the atmosphere and the ocean are calculated herein.
Heat ﬂux The latent and sensible heat ﬂuxes are calculated by using bulk formulae, which
are based on those of Louis (1979) with modiﬁcations by Uno et al. (1995). Roughness ofopen ocean is determined by the method of Miller et al. (1992), and that for snow and ice
covered surface is prescribed at 1.0 × 10−3 m and 2.0 × 10−3 m, respectively. The ﬂuxes
are calculated separately for open water and ice/snow-covered surface at each horizontal
atmosphere-sea ice coupler grid (see section 1.1 for what the “atmosphere-sea ice coupler”
means). The upward longwave radiative ﬂux is calculated by the Stefan-Boltzmann’s fourth
power low of the surface temperature, with the emissivity of 0.95. This calculation is also
made separately for open water and ice/snow-covered surface. The surface albedo of open
water depends on the solar zenith angle. It is calculated by a three-coeﬀcient polynomial ﬁtof the Payne’s (1972) data. Albedo of the bare ice surface is ﬁxed at 0.5. Albedo of the snowcovered surface depends on temperature by taking into consideration the existence of partialsnow-cover at a relatively high temperature. It is 0.9 for temperature lower than 15◦ C, is 0.7
for 0◦ C, and changes linearly in between for HI. For MID, it is 0.85 for temperature lower
than 15◦ C, is 0.65 for 0◦ C, and changes linearly in between.
Together with the downward shortwave and longwave radiative ﬂuxes, which are calculated by the atmospheric model, the air-sea/ice heat ﬂux is determined. For the ice/snow
surface, the surface temperature is also determined in this procedure. Penetration of theshortwave radiative ﬂux into snow or sea ice is not taken into account.
Freshwater ﬂux The freshwater ﬂux due to evaporation or sublimation is determined inconsistent with the latent heat ﬂux. When sea ice sublimates, the contained salt is drained
into the ocean in order to guarantee the salt conservation. Precipitation calculated by theatmospheric model is passed to the sea ice model. Snowfall accumulates on sea ice surface,
while rainfall is directly added to seawater below. When snow falls on open water, thecorresponding latent heat of melt is taken away from the top level. Liquid and solid waterruns oﬀ from the river model. Liquid water runoﬀ is added to rainfall, and solid water runoﬀ
is added to snowfall.

28

Momentum ﬂux The wind stress is also calculated by use of the bulk formula of Louis(1979), with modiﬁcations by Uno et al. (1995) and the open water surface roughness
formualtion of Miller et al. (1992). Surface roughness for snow and ice is 1.0 × 102 m and2.0 × 102 m, respectively. The ice-ocean drag stress, obtained as a result of the sea ice
dynamics calculation, is added to the wind stress for open ocean with the weight of the sea
ice concentration and passed to the ocean.
The wind stress τ is ﬁltered in time when it is passed to sea ice or the ocean. The ﬁlteredvalue τ̃ is calculated by



1 t

τ̃ (t) =τ (t )e(t −t)/T dt ,T −∞where t is time and T is the time-scale of the ﬁltering, which is chosen to be 1 day.

29

(12)

References
Adcroft, A., C. Hill, and J. Marshall, 1997: Representation of topography by shaved cells ina height coordinate ocean model, Mon. Wea. Rev., 125, 2293-2315.
Arakawa, A., and W. H. Schubert, 1974: Interactions of cumulus cloud ensemble with thelarge-scale environment. Part I, J. Atmos. Sci., 31, 671-701.
Arakawa A., and M. J. Suarez, 1983: Vertical diﬀerencing of the primitive equations in sigma
coordinates, Mon. Wea. Rev., 111, 34-45.
Ball, J. T., 1998: An analysis of stomatal conductance, PhD dissertation, Stanford University, 89pp.
Berry, E. X., 1967: Cloud droplet growth by collection. J. Atmos. Sci., 24, 688-701.
Beven, K. J., and M. J. Kirkby, 1979: A physically based, variable contributing area modelof basin hydrology, Hydrol. Sci. Bull., 24, 43-69.
Carpenter, R. L., Jr., K. K. Drogemeier, P. R. Woodward, and C. E. Hane, 1990: Applicationof the piecewise parabolic method to meteorological modelling, Mon. Wea. Rev., 118,
586-612.
Colella, P., and P. R. Woodward, 1984: The Piecewise Parabolic Method (PPM) for GasDynamical Simulations, J. Comput. Phys., 54, 174-201.
d’Almeida, G. A., P. Koepke, and E. P. Shettle, 1991: Atmospheric Aerosols: Global Climatology and Radiative Characteristics, Ed. A. Deepak, 561 pp.
Dukowicz, J. K., and R. D. Smith, 1994: Implicit free-surface method for the Bryan-CoxSemtner ocean model, J. Geophys. Res., 99, 7991-8014.
Emori, S., T. Nozawa, A. Numaguti and I. Uno, 2001: Importance of cumulus parameterization for precipitation simulation over East Asia in June, J. Meteorol. Soc. Japan,
79, 939-947.
Farquhar, G. D., S. von Caimmerer and J. A. Berry, 1980: A biochemical model of photosynthetic CO2 ﬁxation in leaves of C3 species, Planta, 149, 78-90.
Garret, C. J. R., and W. H. Munk, 1975: Space-time scales of internal waves: A progress
report, J. Geophys. Res., 80, 291-297.
30

Geleyn, J. F., and A. Hollingsworth, 1979: An economical analytical method for the computation of the interaction between scattering and line absorption of radiation, Contrib. Atmos. Phys., 52, 1-16.
Gent, P. R., J. Willebrand, T. J. McDougall, and J. C. McWilliams, 1995: Parameterizing eddyinduced tracer transports in ocean circulation models, J. Phys. Oceanogr., 25, 463-474.
Gregg, M. C., T. B. Sanford, and D. P. Winkel, 2003: Reduced mixing from the breaking ofinternal waves in equatorial waters, Nature, 422, 513-515.
Gregory, D., R. Kershaw and P. M. Inness, 1997: Parametrization of momentum transport
by convection. II: Tests in single-column and general circulation models, Q. J. R. Meteorol. Soc., 123, 1153-1183.
Hasumi, H., 2000: CCSR Ocean Component Model (COCO) Version 2.1, CCSR Report, 13,
68pp.
Hibler, W. D., III, 1979: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9,815-846.
Hobbs, P. V., J. S. Reid, R. A. Kotchenruther, R. J. Ferek and R. Weiss, 1997: Direct radiative forcing by smoke from biomass burning, Science, 275, 1776-1778.
Holtslag, A. A. M., and B. Boville, 1993: Local versus nonlocal boundary-layer diﬀusion ina global climate model, J. Climate, 6, 1825-1842.
Hunke, E. C., and J. K. Dukowicz, 1997: An elastic-viscous-plastic model for sea ice dynamics, J. Phys. Oceanogr., 27, 1849-1867.
Ishizaki, H., and T. Motoi, 1999: Reevaluation of the Takano-Oonishi scheme for momentum advection on bottom relief in ocean models, J. Atmos. Oceanic Tech., 16, 19942010.
Killworth, P. D., D. Stainforth, D. J. Webb, and S. M. Paterson, 1991: The development ofa free-surface Bryan-Cox-Semtner ocean model, J. Phys. Oceanogr., 21, 1333-1348.
Leonard, B. P., 1979: A stable and accurate convective modelling procedure based on quadratic
upstream interpolation, Comp. Methods Appl. Mech. Eng., 19, 59-98.
Leonard, B. P., M. K. MacVean, and A. P. Lock, 1993: Positivity-preserving numerical schemes
for multidimensional advection, NASA Tech. Memo. 106055, ICOMP-93-05, 62pp.
31

Leonard, B. P., M. K. MacVean, and A. P. Lock, 1994: The ﬂux-integral method for multidimensional convection and diﬀusion, NASA Tech. Memo. 106679, ICOMP-94-13,
29pp.
Le Treut, H., and Z.-X. Li, 1991: Sensitivity of an atmospheric general circulation model toprescribed SST changes: feedback eﬀects associated with the simulation of cloud opticalproperties, Climate Dynamics, 5, 175-187.
Lin, S.-J., and R. B. Rood, 1996: Multidimensional ﬂux-form semi-Lagrangian transport schemes,Mon. Wea. Rev., 124, 2046-2070.
Lin, S.-J., W. C. Chao, Y. C. Sud and K. Walker, 1994: A class of the van Leer-type transport schemes and its application to the moisture transport in a general circulation
model, Mon. Wea. Rev., 122, 1575-1593.
Louis, J., 1979: A parametric model of vertical eddy ﬂuxes in the atmosphere, Bound.Layer Meteor., 17, 187-202.
Maykut, G. A., 1982: Large-scale heat exchange and ice production in the Central Arctic,
J. Geophys. Res., 87, 7971-7984.
Mellor, G. L., and T. Yamada, 1982: Development of a turbulence closure model for geostrophic
ﬂuid problems, Rev. Geophys., 20, 851-875.
Mellor, G. L., and L. Kantha, 1989: An ice-ocean coupled model, J. Geophys. Res., 94, 10,93710,954.
Mesinger, F., and A. Arakawa, 1976: Numerical methods used in atmospheric models, GARP
Publications Series 17, WMO–ICSU Joint Organising Committee, 64pp.
Miller. M. J., A. C. M. Beljaars, and T. N. Palmer, 1992: The sensitivity of the ECMWF
model to the parameterization of evaporation from the tropical oceans, J. Climate,5, 418-434.
Nakajima, T., M. Tsukamoto, Y. Tsushima, A. Numaguti and T. Kimura, 2000: Modellingof the radiative process in an atmospheric general circulation model, Appl. Opt., 39,4869-4878.
Nakano, H., and N. Suginohara, 2002: Eﬀects of bottom boundary layer parameterization
on reproducing deep and bottom waters in a World Ocean model, J. Phys. Oceanogr.,
32, 1209-1227.
32

Noh, Y., and H.-J. Kim, 1999: Simulations of temperature and turbulence structure of theoceanic boundary layer with the improved near-surface process, J. Geophys. Res., 104,
15,621-15,634.
Numaguti, A., M. Takahashi, T. Nakajima and A. Sumi, 1997: Description of CCSR/NIES
Atmospheric General Circulation Model, CGER’s Supercomputer Monograph Report,
Center for Global Environmental Reserch, National Institute for Environmental Studies, No.3, 1-48.
Oki, T., and Y. C. Sud, 1998: Design of total runoﬀ integrating pathways (TRIP) - a global
river channel network, Earth Interactions, 2.
Pan, D.-M., and D. A. Randall, 1998: A cumulus parameterization with a prognostic closure, Q. J. R. Meteorol. Soc., 124, 949-981.
Payne, R. E., 1972: Albedo of the sea surface, J. Atmos. Sci,, 29, 959-970.

Rosati and Miyakoda, 1988: A general circulation model for upper ocean simulation, J. Phys. Oceanogr.
18, 1601-1626.
Semtner,A. J., Jr., 1976: A model for the thermodynamic growth of sea ice in numericalinvestigations of climate, J. Phys. Oceanogr., 6, 379-389.
Shibata, K., and A. Uchiyama, 1992: Accuracy of the delta-four-stream approximation in
inhomogeneous scattering atmosphere, J. Meteorol. Soc. Japan, 70, 1097-1109.
Smagorinsky, J., 1963: General circulation experiments with the primitive equations, Mon. Wea. Rev.,91, 99-164.
Smith, R. N. B., 1990: A scheme for predicting layer clouds and their water content in a
general circulation model, Q. J. R. Meterol. Soc., 116, 435-460.
Takata, K., T. Watanabe and S. Emori, 2003: Development of the minimal advanced treatments of surface interaction and runoﬀ, Global and Planetary Change, 38, 209-222.
Takemura, T., H. Okamoto, Y. Maruyama, A. Numaguti, A. Higurashi and T. Nakajima, 2000:
Global three-dimensional simulation of aerosol optical thickness distribution of various
origins, J. Geophys. Res., 105, 17,853-17,873.

33

Takemura, T., T. Nakajima, O. Dubovik, B. N. Holben and S. Kinne, 2002: Single-scattering
albedo and radiative forcing of various aerosol species with a global three-dimensionalmodel, J. Climate, 15, 333-352.
Tang, I. N., and H. R. Munkelwitz, 1994: Water activities, densities, and refractive indicesof aqueous sulfates and sodium nitrate droplets of atmospheric importance, J. Geophys. Res., 99, 18,801-18,808.
Tsujino, H., H. Hasumi and N. Suginohara, 2000: Deep Paciﬁc circulation controlled by vertical diﬀusivity at the lower thermocline depths, J. Phys. Oceanogr., 30, 2853-2865.
Uno, I., X.-M. Cai, D. G. Steyn, and S. Emori, 1995 A simple extension of the Louis method
for rough surface layer modelling, Bound.-Layer Meteor., 76, 395-409.
Watanabe, T., 1994: Bulk parameterization for a vegetated surface and its application to asimulation of nocturnal drainage ﬂow, Boundary-Layer Meteorol., 70, 13-35.
WCP-55, 1983: Reports of the experts meeting on aerosols and their climatic eﬀects, Eds. A. Deepakand H. E. Gerber, World Meteorological Organization, 107 pp.
Wiscombe, W. J., and S. G. Warren, 1980: A model for the spectral albedo of snow: I. Puresnow, J. Atmos. Sci., 37, 2712-2733.

34