TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Part V

Boundary Processes

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 14

Sea surface fluxes
The governing equations introduced in Chapter 2 need boundary conditions. At the sea surface, they take the form of theflux boundary condition. Momentum fluxes for the equations of motion, heat and material fluxes for the temperature andother tracer equations, and fresh water fluxes for the continuity equation. These surface fluxes are the main driving forceof oceanic motions. This chapter explains how they are treated in MRI.COM.
Ideally, the realistic momentum, heat and fresh water fluxes based on observation should be used when a realistic
simulation is intended. However, because existing sea surface fluxes derived from observations include large errors andthe ocean models have their own errors, it is not common for an ocean-sea ice model in stand-alone mode to be driven
by fluxes without a feedback process. It is a common practice for an ocean-sea ice model to compute fluxes using a bulk
formula with sea surface meteorological elements derived from atmospheric reanalysis data and the surface temperatureand velocity of its own.
For a simulation using a coupled atmosphere-ocean model, all surface fluxes may be given by the atmospheric component
that computes fluxes using the surface states of the oceanic component. Users are referred to Yukimoto et al. (2011) for
the treatment of fluxes in MRI-ESM.
Section 14.1 describes momentum flux. How heat fluxes are treated is described in Sections 14.2 through 14.5.
Treatment of fresh water fluxes is described in Section 14.6, including the computations of the equivalent fluxes underinvariable first layer volume condition (Section 14.6.2). Sections 14.7 and 14.8 explain how heat and fresh water fluxes
are used in temperature and salinity equations. Several options are available for the choice of the bulk formula to calculate
momentum, latent and sensible heat, and evaporative fluxes. They are detailed in Section 14.9 along with a generalformulation for bulk transfer coeﬃcients. How to give the attributes of external forcing data to the model is explained in
Section 14.10. Finally, some technical details are presented in Sections 14.11 and 14.12.
Here is a cautionary note. Although recent satellite observations enable us to obtain sea surface fluxes with highresolution in space and time, even higher accuracy is necessary for practical uses. For example, a bias of several W m−2
in heat flux may greatly aﬀect the thickness of sea ice, meaning that accuracy on the order of several W m−2 is necessary
to clarify climatic change (WGASF, 2000). Eﬀorts in enhancing observations and evaluating sea surface fluxes based on
various methods have been expanded globally. In the future, a high-resolution ocean model and a new scheme for advection
and diﬀusion may be developed to improve the simulation capability. It is noted, however, that increased observation
frequency does not necessarily guarantee improved accuracy of the fluxes (for example, it is unlikely that the accuracy ofa bulk coeﬃcient would be improved). Hence, it should be kept in mind that sea surface fluxes presently involve large
uncertainties.

14.1 Momentum flux (surface stress)
14.1.1 Formulation
Surface forcing to the momentum equation, or surface momentum flux into the ocean, is in the form of wind stress (orstress from sea ice), and is treated as a body force to the first level velocity in the model algorithm(τx, τy )∂(u1, v1 ) = ... +,∂tρ0 ∆z 1

(14.1)

2

where (u1, v1 ) is horizontal velocity at the first level, ∆z 1 is the thickness of the first layer, ρ0 is the reference density of2
sea water (Table 2.1), and τx and τy are zonal and meridional stresses.
Surface stress is an area weighted average of stresses due to surface wind and sea ice:(τx, τy ) = (1 − A)(τAOx , τAOy ) + A(τIOx , τIOy ),– 149 –

(14.2)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.1 Momentum flux (surface stress)→ and −
where A is the area fraction of sea ice. −τ−AOτ−→IO are air-ocean and ice-ocean stresses, respectively.
−→−−→ is computed by using the wind vector of the surface air, U
Wind stress (air-sea momentum flux) vector, τAO
a , the first−→level velocity (surface current), Us , and a bulk transfer coeﬃcient CD (see Section 14.9 for computation) as follows:
−→
−
→ −→
−
→
−
→ = ρ C |Uτ−AOa D a − αUs |(Ua − αUs ),(14.3)

where ρ a is the density of air and α is the contribution of the surface current to the calculation of relative surface wind.
Usually, α is unity, but it can be tuned to reduce damping eﬀect of the wind stress on the oceanic current (see Table 14.3).
See Chapter 17 for the computation of ice-ocean stress (−
τ−→
IO ).

14.1.2 Numerical implementation
On implementing Eq. (14.3), the time level of the surface velocity must be the starting time level of the integration. For
the time level n of the leap-frog scheme,−→n−→n−1 −→n−→−→ n = ρ C |U
τ−AO|(Ua − αUs n−1 ),a D a − αUs(14.4)

and for the first part of the Matsuno scheme,−→n−→n −→n−→n−→ n = ρ C |Uτ−AOa D a − αUs |(Ua − αUs ).

(14.5)The same treatment should be applied to the computation of ice-ocean stress. The surface velocity at the previous timelevel should be passed to the sea ice part for the leap-frog scheme.
In the bulk formulae implemented in MRI.COM, bulk transfer coeﬃcients are computed on the tracer point. Those on
the velocity point are obtained as area-weighted averages from the surrounding tracer points.

14.1.3 Usagea. Model option

By default, wind stress vector is given as external data and each component should be read from separate files. (Windstress data is always required, and thus files filled with zeros are necessary when wind stress is not applied.) The input data
should be in units of [dyn cm−2 ]. When TAUBULK option is specified, the wind stress is calculated by using a bulk formula.
In this option, components of wind vector at 10 m in units of [cm s−1 ] should be prepared as external data, instead of thedefault wind stress data.
b. Input data

Time series of wind stress or surface wind vector data on a structured lattice should be prepared in the single precision,
direct access format with a constant time interval. The data may be interpolated in time and space in the model. Attributesof the data are specified by nml_force_data. The name of the data and the units assumed by the model are listed on
Table 14.1. The name of the data should be "U-wind" and "V-wind" regardless of wind stress or wind vector. See Section
14.10 for details on how to specify external forcing files.
Table14.1 Wind data to be read from namelist nml_force_data.

X-ward wind stressY-ward wind stress
X-ward wind speed
Y-ward wind speed

nameU-windV-windU-windV-windunits
dyn cm−2
dyn cm−2
cm s−1cm s−1usage
if not TAUBULK
if not TAUBULK
if TAUBULK
if TAUBULK

c. Namelist

At the initial stage of the model integration, wind stress or wind may be gradually set up from zero at the start and unity
after a specified period. The set up period is specified in namelist nml_force_windsetup (Table 14.2).

– 150 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes
Table14.2 Namelist nml_force_windsetup.

variable nameidays_wind_setup

units
day

description
wind set up period in days

usage
default : no tapering on wind

Contribution of the surface current to the calculation of relative surface wind in (14.3) may be tuned by namelist
nml_bulk_wind (Table 14.3). This is because the use of full surface current may lead to too strong damping on swift
surface current in some model settings.
Table14.3 Namelist nml_bulk_wind.

variable namecurrent_contri_factor

units
1

description
Contribution of the surface current to thecalculation of relative surface wind

usage
default = 1.0

14.2 Heat flux
14.2.1 Formulation
There are four main components in surface heat flux: shortwave radiation (QSH ), longwave radiation (QLO ), latent (QLA ),and sensible (QSN ) heat flux. Specific expression for these fluxes is shown in the following sections. In addition to thoseθ ) contribute
fluxes, the heat exchanges with the sea ice (Qice ) and the heat transport due to surface fresh water fluxes (FWF
to surface heat flux. See Section 14.6 for fresh water fluxes and Chapter 17 for details of exchange with sea ice.
The total sea surface heat flux QNET is given byθ
QNET = QSH + QLO + QLA + QSN + Qice + ρ0 cp FWF,(14.6)

where ρ0 and cp are density and specific heat of sea water. Note that downward heat flux is defined as positive.

14.2.2 Numerical implementation
The components of net heat flux (14.6) are broken into three terms,
θ
QNET = QSH + QOTHER + ρ0 cp FWF.

(14.7)

Shortwave radiation (QSH ) penetrates into the ocean interior. Other terms (QOTHER ) are absorbed in the first layer. Freshθ ) needs separate treatment because it should be consistent with the continuity equation andwater transport term (ρ0 cp FWF
included in the solution procedures of the tracer advection term.

14.2.3 Usagea. Model options

The HFLUX option must be selected to treat surface heat fluxes explicitly in the model.
b. Input data

Time series of shortwave and longwave radiation, SST (unnecessary for LWDOWN option), surface air temperature, surface
specific humidity, surface wind speed, and sea-level pressure on a structured lattice should be prepared. The data may beinterpolated in time and space in the model. Attributes of the data are specified by nml_force_data. The name of thedata and the units assumed by the model are listed on Table 14.4. See Section 14.10 for details on how to specify external
forcing files at run time.

– 151 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.3 Shortwave radiation flux
Table14.4 Heat flux data to be read from namelist nml_force_data when HFLUX option is selected.

elements
Downward shortwave radiation
Downward longwave radiation

nameShortWave
LongWave

Sea surface temperature
Surface air temperature
Surface air specific humidity
Scalar wind speed
Sea level pressure

SST
TempAir
SphAirScalarWind
SeaLevelPressureunits
erg s−1 cm−2 = ×10−3 W m−2
erg s−1 cm−2

usage

◦C

Net longwave by default(SST is required).
Downward if LWDOWN
unnecessary if LWDOWN

1
cm s−1hPa

unnecessary if TAUBULK
also available for SLP◦C

14.3 Shortwave radiation flux
14.3.1 Formulation
By default, the downward shortwave radiation (Qdown
SH ) is read as external forcing data. A part of the downward shortwave
radiation (α o Qdown)
isreflected
at
the
sea
surface,and the remainder penetrates into the ocean interior as the net shortwaveSH
downradiation flux (QSH (z = η) = (1 − α o )QSH ), where α o is the albedo at the sea surface and η is the position of sea surface
in the fixed vertical coordinate. More than 50% of the radiation (in the near infra-red band) is absorbed within a depth of
1 m below the sea surface, but the remainder (in the visible and ultraviolet bands) penetrates further into the ocean interior
and aﬀects the subsurface temperature structure.
The penetrating shortwave radiation is expressed as a function of depth (z) in a fixed frame QSH (z) and its verticalconvergence ∂QSH (z)/∂z warms the sea water.
∂θ
1 ∂QSH (z)= ... +,∂t
ρ0 c p∂z(14.8)

where ρ0 is the density of sea water and cp is the specific heat of sea water. There are three options regarding how thevertical profile of shortwave penetration is specified.
a. Standard scheme of penetration

According to Paulson and Simpson (1977), the shortwave radiation flux penetrating into the ocean interior is given byusing the fraction of the near infra-red band (R) and the attenuation distances for the near infra-red band (ζ1 ) and thevisible and ultraviolet band (ζ2 ) as,
QSH (z) = QSH (0)[R exp{(−η + z)/ζ1 } + (1 − R) exp{(−η + z)/ζ2 }],(14.9)

where we set R = 0.58, ζ1 = 35 cm, and ζ2 = 2300 cm, using the optical properties of Water Type I based on theclassification by Jerlov (1976). QSH (0) is the net solar radiation flux at the sea surface.
b. A penetration scheme that takes into consideration of the insolation angle

The SOLARANGLE option enables us to include the eﬀect of the insolation angle on the shortwave radiation flux in theoceanic interior. This scheme is based on Ishizaki and Yamanaka (2010). In this case, if the depth (z) on the r.h.s. of
Equation (14.9) is replaced by the penetrating distance from the sea surface, the shortwave radiation is given by
QSH (z) = QSH (0)[R exp{(−η + z)/(ζ1 sin θ w )} + (1 − R) exp{(−η + z)/(ζ2 sin θ w )}],

(14.10)where θ w is the penetrating angle in the ocean interior. When SOLARANGLE option is specified in coupled models, Equation(14.10) is replaced by
QSH (z) = QSHb (z) + QSHd (z),(14.11)
where QSHb (z) and QSHd (z) are the shortwave radiation fluxes due to direct and diﬀuse solar radiation. Those fluxes are
expressed as follows:
QSHb (z) = QSHb (0)[Rb exp{(−η + z)/(ζ1 sin θ w )} + (1 − Rb ) exp{(−η + z)/(ζ2 sin θ w )}],
QSHd (z) = QSHd (0)[Rd exp{(−η + z)/ζ1 } + (1 − Rd ) exp{(−η + z)/ζ2 }],– 152 –

(14.12)(14.13)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes

where Rb and Rd are the ratios of near infra-red to total radiation in direct and diﬀuse solar radiation, respectively. QSHb (0)and QSHd (0) are the net direct and diﬀuse solar radiation fluxes at the sea surface.
c. A penetration scheme that takes into consideration of chlorophyll concentration

Recent studies indicate that solar radiation absorption and local heating within the upper ocean are strongly influenced by
the chlorophyll concentration. By specifying CHLMA94 and NPZD options, we can use a shortwave penetration model withthe chlorophyll concentration (Morel and Antoine, 1994). In this scheme, the shortwave radiation flux penetrating into theocean interior is given by
QSH (z) = QSH (0)[R exp{(−η + z)/ζ0 } + (1 − R)(V1 exp{(−η + z)/ζ1 } + V2 exp{(−η + z)/ζ2 })].

(14.14)

The first exponential is for the infra-red waveband (> 750 nm), which is not influenced by biological materials andattenuates with a distance of ζ0 = 26.7 cm. The second and third exponentials are for the visible and ultraviolet bands (<750 nm). V1 , V2 , ζ1 , and ζ2 are calculated from an empirical relationship as a function of chlorophyll concentration (chl[mg · m−3 ]) as follows:
V1 = 0.321 + 0.008C + 0.132C 2 + 0.038C 3 − 0.017C 4 − 0.007C 5,(14.15)

V2 = 0.679 − 0.008C − 0.132C 2 − 0.038C 3 + 0.017C 4 + 0.007C 5,(14.16)5

ζ1 = 1.540 − 0.197C + 0.166C − 0.252C − 0.055C + 0.042C ,(14.17)

ζ2 = 7.925 − 6.644C + 3.662C − 1.815C − 0.218C + 0.502C ,(14.18)

2

3

2

3

4

4

5

where C = log10 (chl). It is noted that V1 + V2 = 1.
When SOLARANGLE option is added to the options mentioned above, the shortwave radiation is slightly modified by
QSH (z) = QSH (0)[R exp{(−η + z)/(ζ0 sin θ w )}+ (1 − R)(V1 exp{(−η + z)/(ζ1 sin θ w )} + V2 exp{(−η + z)/(ζ2 sin θ w )})].

(14.19)In coupled models, Equation (14.19) is replaced by
QSH (z) = QSHb (z) + QSHd (z),(14.20)where QSHb (z) and QSHd (z) are the shortwave radiation fluxes due to direct and diﬀuse solar radiation. These fluxes are
expressed as follows:
QSHb (z) = QSHb (0)[Rb exp{(−η + z)/(ζ0 sin θ w )}+ (1 − Rb )(V1 exp{(−η + z)/(ζ1 sin θ w )} + V2 exp{(−η + z)/(ζ2 sin θ w )})],
QSHd (z) = QSHd (0)[Rd exp{(−η + z)/ζ0 }+ (1 − Rd )(V1 exp{(−η + z)/ζ1 } + V2 exp{(−η + z)/ζ2 })],

(14.21)(14.22)where Rb and Rd are the ratios of near infra-red to total radiation for direct and diﬀuse solar radiation, respectively.
d. Albedo schemes

There are three options for the sea surface albedo. The first option (albedo_choice = 1 in namelist nml_albedo_ocean)
is a constant value, which should be specified as alb in namelist nml_albedo_ocean. The second option(albedo_choice = 2) is based on Large and Yeager (2009) and is given by
α o = 0.069 − 0.011 cos(2ϕ),(14.23)where ϕ is latitude. The third option (albedo_choice = 3) is based on Baker and Li (1995) and is given by
α o = 0.06 + 0.0421x 2 + 0.128x 3 − 0.04x 4 +(0.074x ) 5
3.12
+
x ,
5.68 + U 1.0 + 3.0U

where x = 1 − sin θ a (θ a is a height angle of the Sun), and U is the surface wind speed [m s−1 ].

– 153 –

(14.24)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.3 Shortwave radiation flux

14.3.2 Numerical implementation
Shortwave flux is computed on T-cells, which means that variations of horizontal cross section due to topography (S) must
be taken into account. The vertical profile of the shortwave penetration can be expressed as a function of the vertical depth
from the sea surface. Discrete form of the absorption of shortwave heat for a grid cell (i, j, k − 12 ) is
(
)
n
n−1
n
(θV )i,n+1
1 = (θV )1 + ... + 2∆t Q SH i, j,k−1 Si, j,k− 1 − Q SH i, j,k Si, j,k+ 1 ,j,k−i, j,k−2

2

2

2(14.25)for the leap-frog scheme.

14.3.3 Usagea. Model option• SOLARANGLE must be selected to take the solar penetration angle into consideration• CHLMA94 must be selected to take the eﬀect of chlorophyll concentration into consideration. The chlorophyll
concentration will be taken from either the biogeochemical model with NPZD option or climatology (name =ChlorophyllClimatology) with CHLCLIM option.
When CHLCLIM option is selected, time series of 2-D chlorophyll concentration should be prepared. The data may beinterpolated in time and space in the model. See Section 14.10 for how to specify external forcing files.
b. Namelist

Namelists related to shortwave radiation are listed on Tables 14.5 through 14.7.
Table14.5 Namelist nml_albedo_ocean.

variable name
alb

units
no (real8)

albedo_choiceno (integer)

description
constant albedo
for albedo_choice = 1
choice of albedo schemeusage

1 : constant given by alb
2 : Large and Yeager (2009)
3 : MRI-ESM1

Table14.6 Namelist nml_solarangle.

variable nameinterval_shortwave_sec

units
sec

l_variable_zenith_ang

logicall_const_zenith_ang_mean

logicall_const_zenith_ang_noon
l_diurnal_var_shortwave

logicallogicall_penetr_ang_considered

logicaldescription
time interval of shortwave data
use zenith angle withdaily variation
use average zenith angle
from sunrise to sunset
use zenith angle at noon
daily variational amplitude of solar flux is calculated in this module
introduction of transmitted angle based on
Snell’s law

usage
requireddefault = .true.
valid if l_variable_zenith_ang = .false.
valid if l_variable_zenith_ang = .false.
default = .false., If this is true,
interval_shortwave_sec = 86400
default = .true.

Table14.7 Namelist nml_frac_nearir.

variable nameqsratio_nearIR

units
real(8)description
ratio of near infra-red band to the total
shortwave flux (R)– 154 –

usage
default = 0.58

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes

14.4 Longwave radiation flux
14.4.1 Formulation
Net longwave radiation into the ocean is the residual between the downward longwave radiation Qdown
LO (external data) andthe upward black body radiation of sea water with the sea surface temperature (θ surf ),
down
4
Qnet
LO = Q LO − em σ(θ surf + 273.15) .

(14.26)

Here, em is the emissivity for sea water, which should be specified by user, and σ = 5.67 × 10−5 erg · s−1 · cm−2 · K−4 isthe Stefan-Boltzmann constant.

14.4.2 Numerical implementation
There is no special issue on the numerical implementation. The net longwave flux is absorbed in the first layer of themodel. Before used to update the first level temperature, it is summed into an array (qothr) which gathers all heat flux
components except for shortwave. Note that the sea surface (first model level) temperature at the starting time level is used
to evaluate upward longwave flux in any time integration scheme.

14.4.3 Usagea. Model option

By default, the net longwave radiation is read as external forcing data. By choosing LWDOWN option, the downward
longwave radiation is read. LWDOWN option will become a standard (default) choice in the near future.
b. Namelist

Emissivity (em ) of black body radiation from the sea surface may be specified by user at run time with namelist
nml_emissivity_sea_water (Table 14.8).
Table14.8 Namelist nml_emissivity_sea_water.

variable nameemissivity

units
no (real8)

description
emissivity of upward longwave radiationfrom sea water

usage
default = 1.0

14.5 Latent and sensible heat fluxes
14.5.1 Formulation
The bulk method is used to calculate latent and sensible heat fluxes. In the bulk method, latent heat flux (QLA ) and sensible
heat flux (QSN ) are calculated using bulk transfer coeﬃcients CE and CH (see Section 14.9),⃗ |(qs − qa ),QLA = − ρ a LCE |∆U(14.27)⃗ |(Ts − θ a ).
QSN = − ρ a cpa CH |∆U(14.28)Here, ρ a is the air density, L is the latent heat of evaporation, qa is the specific humidity of the surface air, qs is thesaturated specific humidity of the sea surface temperature, θ a is the surface air potential temperature (K), Ts is the absolute⃗ = U⃗a − αU⃗s is the vector diﬀerence between the surface wind and the first levelsea surface water temperature (K), and ∆Uvelocity (surface current). As explained in relation to the wind stress (Section 14.1), α is the contribution of the surface
current to the calculation of the relative wind vector. The quantity cpa is the specific heat of air.

14.5.2 Numerical implementation
There is no special issue on the numerical implementation. The latent and sensible heat fluxes are absorbed in the firstlayer of the model. Before used to update the first level temperature, they are summed into an array (qothr) which gathers
all heat flux components except for shortwave. Note that the sea surface (first model level) temperature at the starting timelevel is used to evaluate fluxes in any time integration scheme.
– 155 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.6 Freshwater flux

14.5.3 Usage
Model option

MRI.COM can provide two calculation methods for bulk transfer coeﬃcients. One of them must be specified when HFLUX
option is chosen.
• BULKNCAR for Large and Yeager (2004; 2009)• BULKKONDO2 for Kondo (1975)
See Sections 14.9.2 and 14.9.3 for details.

14.6 Freshwater flux
14.6.1 Formulation
Freshwater flux through the sea surface (F W ) is caused mainly by precipitation (P), evaporation (E), river discharge (R),and formation-melting of sea ice (I):F W = P − E + R + I,
(14.29)Precipitation and river discharge are given as input data. Fresh water flux due to sea ice formation-melting is given by thesea ice model. Evaporation is computed by a bulk formula in the ocean model and it should be consistent with the latent
heat flux:⃗ |(qs − qa ) = −QLA /L.
E = ρ a CE |∆U(14.30)Note that evaporation (E) is defined as positive upward.
Fresh water flux is treated as the surface boundary condition for the surface height equation:∂η= ... + F W .
∂t(14.31)

Fresh water is added to/subtracted from the surface grid cell. Thickness of grid cells in the whole vertical column should
be modified to adapt to this change in the surface height. Concentration of materials (tracer species) in a grid cell should
be updated by taking into account of these volume changes. Knowledge about the contents of materials in the fresh water
flux is required for the computation, which is explained below.

14.6.2 Material fluxes due to fresh water flux
Precipitation, evaporation, and river discharge (WAO = P − E + R) are assumed to have the first level temperature (θ surf )and zero salinity. The water transported through the ice bottom (WIbot ) is also assumed to have the first level temperature.
The water transported from the ice surface (WIsurf > 0) is assumed to have the freezing point temperature (θ freeze ). The
water exchanged with ice is assumed to have low salinity (SI = 4.0 pss). Note that the freezing point temperature (θ freeze )is given by mSI , where m = −0.0543 ◦C/pss.
a. When fresh water fluxes are explicitly added to the systemIf fresh water flux is simply added to the first layer volume, surface temperature and salinity fluxes due to the fresh water
flux are given as follows:θ
FWF= (WAO + WIbot ) · θ surf + WIsurf · θ freeze,(14.32)S
FWF= (WIbot + WIsurf ) · SI .

(14.33)b. When fresh water fluxes are not explicitly added to the system

If the surface fresh water flux is not added to the first layer to avoid a long-term drift of sea level (SFLUXW option),equivalent temperature and salinity fluxes should be imposed to incorporate the eﬀects of surface fresh water flux. To
obtain a formula, the conservation equations for the first layer heat and salinity are considered. Temperature and salinity
are assumed to evolve from the old values (θ old, Sold ) to the new values (θ new, Snew ) due to fresh water flux.

– 156 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes

If a volume change is allowed,
Vnew θ new = Vold θ old + (WAO θ old + WIbot θ old + WIsurf θ freeze ) · ∆A · ∆t,Vnew Snew = Vold Sold + (WIbot SI + WIsurf SI ) · ∆A · ∆t,Vnew = Vold + (WAO + WIbot + WIsurf ) · ∆A · ∆t,(14.34)(14.35)(14.36)

where ∆A is the horizontal area.
If a volume change is not allowed,
θVold θ new = Vold θ old + FWF· ∆A · ∆t,(14.37)

S
Vold Snew = Vold Sold + FWF· ∆A · ∆t,(14.38)θ and F S are temperature and salinity flux for the constant volume case.
where FWF
WF
Removing θ new , Snew , and Vnew from the above equations, we have,

Vold {Vold θ old + (WAO θ old + WIbot θ old + WIsurf θ freeze ) · ∆A · ∆t}θ= {Vold + (WAO + WIbot + WIsurf ) · ∆A · ∆t}(Vold θ old + FWF· ∆A · ∆t),
Vold {Vold Sold + (WIbot SI + WIsurf SI ) · ∆A · ∆t}
S= {Vold + (WAO + WIbot + WIsurf ) · ∆A · ∆t}(Vold Sold + FWF· ∆A · ∆t).

(14.39)(14.40)The fluxes for the constant volume case are given by
Vold
WIsurf · (θ old − θ freeze ),Vold + W · ∆A · ∆t
Vold
S
FWF
=−
{WAO · Sold + (WIbot + WIsurf ) · (Sold − SI )},Vold + W · ∆A · ∆t

θ
FWF
=−

(14.41)(14.42)where W = WAO + WIbot + WIsurf .
If |Vold | ≫ |W · ∆A · ∆t|, we have,θFWF= − WIsurf · (θ old − θ freeze ),(14.43)S
FWF= − WAO · Sold − (WIbot + WIsurf ) · (Sold − SI ).

(14.44)14.6.3 Numerical implementation
a. Adding/subtracting new sea water

Fresh water flux is applied to the sea level prediction in the baroclinic time interval:
η n+1 = η n−1 + ... + 2∆t(F W ) n + (F W ) n−1,
2(14.45)where (F W ) n−1 and (F W ) n are fresh water fluxes evaluated at time levels n − 1 and n, respectively.
In the same manner, fresh water transport is reflected in the equations for temperature(θV ) n+1
= (θV ) n−1+ ... + ∆t [{(WAO ) n−1 + (WIbot ) n−1 } · (θ surf ) n−1 + (WIsurf ) n−1 · (θ freeze ) n−1 ]
1
1
2

2

+ ∆t [{(WAO ) n + (WIbot ) n } · (θ surf ) n + (WIsurf ) n · (θ freeze ) n ](14.46)and salinity(SV ) n+1
= (SV ) n−1+ ... + ∆t {(WIbot ) n−1 + (WIsurf ) n−1 } · SI + ∆t {(WIbot ) n + (WIsurf ) n } · SI .
1
1
2

2

– 157 –

(14.47)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.6 Freshwater flux

b. Suppressing long-term drift of sea water volume

Because the global mean of F W should not necessarily be zero and there is no feedback mechanism in an ocean-sea ice
system to balance the fresh water budget, the globally averaged volume of the ocean may continue to increase or decrease.
To avoid this, the integral of the globally averaged freshwater flux is set to zero by selecting WADJ option,W
FADJ= FW − FW ,(14.48)where F W is the global mean of F W (= P − E + R), where E includes sublimation over sea ice. Note that exchange withsea ice is not included in the computation of F W .
Another option for suppressing the long term drift of sea level is to transform all freshwater fluxes to salinity flux, whichis activated by specifying SFLUXW option. Formulae for the temperature and salinity fluxes with this option are given by(14.43) and (14.44).

14.6.4 Usagea. Model options• WFLUX should be specified to explicitly treat fresh water fluxes in the model. In this case, precipitation data must
be prepared as input data.
• RUNOFF enables river discharge to be treated. River run-oﬀ data should be given as external data.
• WADJ imposes conservation of sea water volume of the ocean - sea ice system.
• SFLUXW computes equivalent temperature and salinity fluxes due to fresh water flux
b. Input data

When WFLUX and RUNOFF options are specified, time series of precipitation and river run-oﬀ data must be prepared,
respectively. Data on a structured lattice should be prepared in the single precision direct access format with a constanttime interval. Attributes of the data are specified by nml_force_data. Name and units assumed by the model aresummarized in Table 14.9.
The data may be interpolated in time and space in the model. See Section 14.10 for how to specify external forcing files.
Table14.9 Water flux data to be read from namelist nml_force_data.

element
precipitation
river discharge

namePrecipitation
RiverDischargeRateunits
g s−1 cm−2
g s−1 cm−2

usage
WFLUX
RUNOFFc. Namelist

The river run-oﬀ data around Antarctica may be treated as the discharge of iceberg if sea ice is present in the grid. Namelist
nml_mkflux specifies this behavior (Table 14.10).
Table14.10 Namelist nml_mkflux. See also Table 14.12.

variable namel_glacierunits
logicall_after_rivermouth_spreadlogicaldescription
if a rivermouth grid has sea ice then riverdischarge is converted into iceberg input
from glacierriver discharge is converted into iceberg after a smoothing operation of river discharge
data

usage
default = .false.
default = .false.

When WADJ option is selected, the global summation of surface fresh water fluxes, which is used for determining theoﬀsetting flux, can be made MPI parallel independent by choosing l_strict_waterflux_adjust = .true. in namelist
nml_waterflux_adjust (Table 14.11).

– 158 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes
Table14.11 Namelist nml_waterflux_adjust.

variable namel_strict_waterflux_adjust

units
logicaldescription
water flux adjustment is made MPI parallel
independent

usage
default = .false.

14.7 Treatment in the temperature equation
14.7.1 Formulation
The contribution of surface forcing (Fzθ ) and shortwave radiation flux (QSH ) to the first layer temperature is expressed asFzθ∂θ
1 ∂QSH
+,
1 = ... +∂t k= 2∆z 1
ρ0 cp ∂z(14.49)2

where ∆z 1 is the thickness of the first layer of a T-cell. The term (Fzθ ) consists of the net sea surface heat flux except2θ ), and restoring of SST (θfor the shortwave radiation flux (QOTHER ), heat transport due to fresh water flux (FWF
surf ) to a∗
specified value (θ ):QOTHER1
θFzθ =(14.50)+ FWF−(θ surf − θ ∗ )∆z 1 .
2
ρ0 c pγθ
The heat flux (QOTHER ) has been converted to temperature transport using the reference density (ρ0 ) and specific heat (cp )of sea water. The third term on the r.h.s. is a restoring term to a specified SST (θ ∗ ). The parameter γθ is a restoring time(in units of seconds) and should be read from input file (see Chapter 13). The restoring term is converted to temperaturetransport by multiplying the thickness of the first layer of a T-cell (∆z 1 ).
2
By default, surface temperature forcing consists only of the restoring term. By specifying HFLUX option, sea surface
heat fluxes are treated explicitly. The heat transport due to fresh water flux are included when WFLUX option is specified.

14.7.2 Numerical Implementation
The heat transport due to fresh water flux is incorporated in the advection schemes. Shortwave absorption is treatedseparately from other heat flux terms. Heat flux components except for the shortwave radiation are gathered in an array(qothr) and added to the tendency of the first layer temperature.
Treatment of supercooled water needs some explanation. When the temperature of the first layer is below the freezing
point, the temperature is set to the freezing point and this heat comes from the latent heat release of the new sea ice
formation. This is done in the sea ice model. Because the sea ice part uses the forward scheme to proceed in time, thesea surface temperature at the "current" time level must be given to the sea ice part from the ocean part. The result is thatthe temperature below the freezing point at the "current" time level is modified and adjusted to the freezing point by thesea ice part. This treatment is not compatible with the leap-frog time stepping employed by the ocean part. Specifically,treating this temperature adjustment process in a form of heat flux exchange may result in a numerical instability. Though
the modification of sea surface temperature in the sea ice part is inconvenient for some purposes, we keep this treatment
until a revision of time-stepping scheme is considered.

14.8 Treatment in the salinity equation
14.8.1 Formulation
The fresh water flux (FW ) modifies the volume but not the salt content of the surface layer, changing the salinity of thesurface layer (see also Chapter 6). Since fresh water flux is not related directly to the sea surface salinity, the model sea
surface salinity might be far from the observed value. Hence, an adjustment is sometimes needed to restore the model sea
surface salinity to the observed one (see the last term on the r.h.s. of the next equation). Generally, the model salinity is
restored to the observed climatological sea surface salinity since no reliable data set of historical sea surface salinity isavailable at present.
The explicit surface forcing (FzS ) consists of salt transport due to formation and melting of sea ice and restoring of SSS(Ssurf ) to a specified value (S ∗ ):FzS = (WIbot + WIsurf ) · SI −– 159 –

1
(Ssurf − S ∗ )∆z 1 .
2
γs(14.51)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.8 Treatment in the salinity equation
The water exchanged with ice (WIbot + WIsurf ) is assumed to have low salinity (SI = 4.0 pss). The parameter γs is a
restoring time (in units of seconds) and should be read from input file (see Chapter 13). The restoring term is converted
to salinity transport by multiplying the thickness of the first layer of a T-cell (∆z 1 ).
2
By default, surface salinity forcing consists only of the restoring term. When WFLUX option is specified, the eﬀect of
the fresh water flux on salinity may be included.
With WFLUX option, by default, the restoring term is converted to the fresh water flux. The expression for F W becomes
FW = P − E + R + I +1 Ssurf − S ∗∆z 1 .
2γs Ssurf(14.52)In this case, the salinity flux should be
FzS = (WIbot + WIsurf ) · SI .

(14.53)

However, the restoring term is not converted to the freshwater flux if SFLUXR option is selected.
If SFLUXW option is specified, the fresh water flux is not added to the ocean and its eﬀect on salinity is now explicitly
expressed as the salinity flux:
FzS = −WAO · Ssurf − (WIbot + WIsurf ) · (Ssurf − SI ) −1(Ssurf − S ∗ )∆z 1 ,2
γs(14.54)

where Ssurf is the first level salinity.
The contribution of surface forcing (FzS ) to the first layer salinity is expressed as
FzS
∂S
,
=
...
+
∂t k= 21∆z 1

(14.55)2where ∆z 1 is the thickness of the first layer of a T-cell.
2

14.8.2 Numerical Implementation
The salinity transport due to fresh water flux is incorporated in the advection schemes. The first layer salinity is calculated
by taking into consideration of the volume change of the grid cell due to the continuity equation.

14.8.3 Usagea. Model option• WFLUX treats fresh water fluxes explicitly in the model. Salinity will be aﬀected by fresh water fluxes.
• SFLUXW converts fresh water flux to surface salinity flux
• SFLUXR suppresses conversion of salinity restoring to water flux
• STABLERUN restores sea surface salinity to climatology when salinity is lower than 5 pss
b. Namelist

Namelists related to salinity flux at the sea surface are listed on Tables 14.12 and 14.13.
Table14.12 Namelist nml_mkflux. See also Table 14.10.

variable namers_localunits

description
how salinity for iceberg and snow is extractedusage
1 : taken from the local grid0 : taken from global sea surface

Table14.13 Namelist nml_sss_restore.

variable namesss_rst_diff_max

units
pss

description
the maximum of the diﬀerence between modeland reference sea surface salinity

usage
default = 100 pss
Continued on next page

– 160 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes

variable namel_nosss_rst_seaicel_nosss_rst_csticel_sss_rst_cnsv
l_strict_sss_rst_cnsv

Table 14.13 – continued from previous page
description
no surface salinity restoring below sea ice
no surface salinity restoring at the coastal gridwith sea ice
logical set globally integrated surface salinity restoringto zerological global summation of surface salinity restoring
flux is done serially so that the result becomes
MPI parallel independentunits
logicallogicalusage
default = .false.
default = .false.
default = .false.
default = .false.

14.9 Bulk transfer coeﬃcientThis section briefly describes how to calculate sea surface fluxes using a bulk formula. For details, readers are referred toKantha and Clayson (2000) and Large and Yeager (2004). Bulk transfer coeﬃcient for air-ice interface is given in Section
17.2.3.

14.9.1 Formulation of bulk formula
Transfer processes through atmosphere and ocean boundaries are governed by complicated turbulent processes. Traditionally, these turbulent fluxes are parameterized as a bulk transfer law. This method does not explicitly solve the interaction
between the atmosphere and the ocean, and attributes all unknown processes to bulk transfer coeﬃcients.
Momentum (τ), sensible heat (HS ), latent heat (H L ), and water vapor (E) fluxes are written in terms of turbulent
components as follows:τ = − ρ a wu = ρ a u∗2,

(14.56)HS = − ρ a cpa wθ = ρ a cpa u∗ θ ∗,(14.57)

H L = − ρ a L E wq = ρ a L E u∗ q∗,(14.58)E = − ρ a wq = ρ a u∗ q∗ = H L /L E ,(14.59)where a bar over two variables denotes a covariance between the turbulent component of vertical velocity (w) and that ofeach physical property, wind speed (u), air potential temperature (θ), and specific humidity (q). ρ a is air density, u∗ isfriction velocity, θ ∗ ≡ HS /( ρ a cpa u∗ ) is the temperature scale in the boundary layer, q∗ ≡ H L /( ρ a L E u∗ ) = E/( ρ a u∗ ) is
the scale of specific humidity in the boundary layer, cpa is the specific heat of air, and L E is the latent heat of water vapor
evaporation.
Sea surface fluxes are also represented using bulk formulae as follows:τ = ρ a CDh |U⃗a − U⃗s | 2,

(14.60)HS = − ρ a cpa |U⃗a − U⃗s |CHh (Ts − θ a ),(14.61)

H L = − ρ a L E |U⃗a − U⃗s |CEh (qs − qa ),(14.62)E = ρ a |U⃗a − U⃗s |CEh (qs − qa ),(14.63)where U⃗a is the wind vector at a height of zh , U⃗s is the current vector at the sea surface, and Ts is the absolute sea surface
water temperature (K). The subscript "a" means the value at a height of z = zh , and the subscript "s" means the value atthe sea surface. The contribution factor (α) of the surface current to the calculation of the relative wind vector explainedin Section 14.1 is omitted here for brevity.
Parameter CDh is called a drag coeﬃcient. CHh and CEh are transfer coeﬃcients for heat and water vapor and are called
the Stanton coeﬃcient and the Dalton coeﬃcient, respectively. These can be estimated by observed atmospheric elements(wind velocity, air temperature, and humidity) at a height in the boundary layer, not by observed turbulent fluxes, using– 161 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.9 Bulk transfer coeﬃcientthe following Equations (14.64), (14.65), and (14.66) and the similarity law of Monin-Obukhov mentioned below,
CDh =

u∗2

,
|U⃗a − U⃗s | 2
HS /ρ a cpa
u∗ θ ∗
CHh = −=,|U⃗a − U⃗s |(Ts − θ a )|U⃗a − U⃗s |(Ts − θ a )
u∗ q∗
E/ρ a
=
.
CEh =
|U⃗a − U⃗s |(qs − qa )|U⃗a − U⃗s |(qs − qa )(14.64)(14.65)(14.66)The similarity law of Monin-Obukhov assumes that physical properties in the atmosphere-ocean boundary layer (a layer
with a thickness of several tens of meters located below the lower mixed layer of the atmosphere) have similar verticalprofiles when they are scaled with the stability parameter and sea surface fluxes. Vertical profiles of wind speed (U), airtemperature (θ), and humidity (q) can be written as follows:
(z)
κz ∂U
= ϕM
,
u∗ ∂z
L
(z)
κz ∂θ
= ϕH
,
θ ∗ ∂z
L
(z)
κz ∂q
= ϕE
,
q∗ ∂z
L

(14.67)(14.68)(14.69)where κ = 0.4 is the von Karman constant, and L is the Monin-Obukhov length scaleL=u∗3u3 θ v.
=− ∗κQ b
κ g wθ v(14.70)In (14.70), Q b = −g wθ v /θ v is the buoyancy flux, θ v is the virtual temperature (θ v = θ (1 + 0.6078 q)), and wθ v is thenet heat flux including the water vapor flux:wθ v = wθ (1 + 0.6078 q) + 0.6078 θ wq.

(14.71)

In (14.67) to (14.69), ζ = z/L is the Monin-Obukhov similarity variable, and ϕ M, H,E is a nondimensional function forwind velocity, air temperature, and specific humidity. The nondimensional function is assumed to be a mathematicallysimple function.
Integrating Equations (14.67), (14.68), and (14.69) vertically, we have the following.
[]u∗z
U (z) − Us =ln − ΨM (ζ ) ,(14.72)
κ
z0
[
]
θ∗
z
θ(z) − Ts =ln
− ΨH (ζ ) ,(14.73)
κ
z0T
[
]
z
q∗
q(z) − qs =ln
− ΨE (ζ ) ,(14.74)κz0E
where
ΨM, H,E (ζ ) =

∫ ζ[

′

1 − ϕ M, H,E (ζ )0] dζ ′
ζ′.

(14.75)

In (14.72) to (14.74), z0 , z0T , and z0E are roughness lengths for each physical property. When the stability of the boundary
layer is already known, bulk transfer coeﬃcients can be estimated using these roughness lengths,
CDh = [

CHh,Eh = [
ln

κ2
zh
ln − ΨM (ζh )z0zh− ΨM (ζh )z0][– 162 –

]2,

κ2
ln

zh
z0T,0E(14.76)].
− ΨH,E (ζh )(14.77)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes

In neutral stability (ζ = 0), the bulk transfer coeﬃcient is a function of the roughness length (z0N ,z0TN , and z0EN foreach) only and is expressed as follows:κ2
CDNh = [(14.78)
]2,
zh
ln
z0N

CHNh,ENh = [

zh
ln
z0N

][

κ2
zh

ln

]

(14.79)z0TN,0EN
CDNh]=[
z0TN,0EN
1 1/2
1 − CDNh ln
κ
z0N
1/2
κCDNh

=[
ln

(14.80)].

zh(14.81)z0TN,0EN

Normally, these neutral bulk transfer coeﬃcients are estimated at a height of 10 m. Non-neutral bulk transfer coeﬃcientsat an arbitrary height (zh ) are connected with the neutral bulk transfer coeﬃcients at a height of 10 m by eliminating theroughness length as follows:
CDN10
CDh = [
(14.82)()] 2 ,
1zh−1/2
1 + CDN10 CDN10 ln− ΨM (ζh )
κ
z10
) 1/2
CDh
CDN10
(
)] ,
CHh,Eh = [
1
zh
−1/2
ln
− ΨH,E (ζh )1 + CHN10,EN10 CDN10κz10(

CHN10,EN10(14.83)where z10 means z = 10 m.
The neutral bulk transfer coeﬃcients at a height of 10 m (CDN10 , CHN10 , and CEN10 ) are often estimated, according to
the stability, as a function of velocity at 10 m. Various formulations can be used (see the following subsections). To be
more realistic, the bulk transfer coeﬃcients should be regarded as a function of wave age, but its general formulation has
not been achieved yet.
In principle, if we have a complete set of atmospheric elements in the boundary layer, we can get the bulk coeﬃcientsand fluxes because height can be adjusted using the similarity law. This requires an iteration of the calculation as follows(you can select an appropriate way for your data):• Convert wind speed, air temperature, and specific humidity into those at 10 m in the neutral stability and use a bulk
coeﬃcient at 10 m in the neutral stability.
• Convert a bulk coeﬃcient at 10 m into one with a height and a stability at which atmospheric elements are observed.
• Convert air temperature and specific humidity into those at a height where wind speeds are observed, and convert
the bulk coeﬃcient at 10 m in the neutral stability into one with a height and a stability at which wind speeds areobserved.

14.9.2 Large and Yeager (2004; 2009): BULKNCAR
The bulk formula based on Large and Yeager (2004; 2009) closely follows the formulation explained in the previous
subsection. Each flux is estimated at a height where wind speed is observed by transforming air temperature and humidityfrom the height where they are observed (z θ and zq , respectively) to the height where wind speed is observed. To do thisiteration in the program bulk.F90, specify BULKITER option in configure.in. The observed heights are specified in
namelist nml_srfflx_dbase (see Table 14.18).
The bulk transfer coeﬃcient at 10 m in the neutral stability is given as follows:{
10 CDN10 =
3

U10N
2.70
−10U 6 ,U10N + 0.142 + 13.09 − 3.14807 × 10
10N

2.34
– 163 –

U10N < 33 m s−1
U10N ≥ 33 m s−1(14.84)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.9 Bulk transfer coeﬃcient√
103 CHN10,EN10 = 18.0 CDN10, stable ζ > 0,√
103 CHN10,EN10 = 32.7 CDN10, unstable ζ ≤ 0,

(14.85)(14.86)

Each physical property is estimated at a height where wind speed is observed in the following operations. First, calculate
the virtual temperature θ v as follows:θ v = θ(z θ )(1 + 0.6078 q(zq ))(14.87)
Next, calculate the bulk transfer coeﬃcient at 10 m in the neutral stability assuming that the first guess for the 10 m wind
⃗ (zu )| = |U⃗a (zu ) − U⃗s | (Equations (14.84) to (14.86)). The first guesses of the
speed in the neutral stability is U10N = |∆U
scales for the friction velocity, air temperature, and specific humidity are estimated assuming that these bulk coeﬃcientsare at the observed height and stability,√√τ⃗ (zu )|,
= CDN10 |∆U
ρa
HS
CHN10
θ∗ =
=√
(θ(z θ ) − Ts ),ρ a c p u∗
CDN10
E
CEN10
q∗ =
(q(zq ) − qsat (Ts )),=√ρ a u∗
CDN10u∗ =(14.88)(14.89)(14.90)where qsat (Ts ) is the saturated specific humidity at sea surface temperature Ts .
Next, perform the iteration using the three Monin-Obukhov similarity variables, ζu = zu /L, ζ θ = z θ /L, and ζ q = zq /L,and an integral of the nondimensional profile function for the vertical gradient of each physical property, ΨM (ζ ) for
momentum and ΨH (ζ ) for scalars.
The Monin-Obukhov similarity variables are calculated as follows:][κgz θ ∗q∗
ζ= 2+.
(14.91)u∗ θ v (q(zq ) + 0.6078−1 )The integral of the non-dimensional profile function is expressed asΨM (ζ ) = ΨH (ζ ) = −5ζ,

(14.92)

if it is stable (ζ ≥ 0), and
(

)
(
)
1+X
1 + X2πΨM (ζ ) = 2 ln
+ ln
− 2 tan−1 (X ) + ,
2
2
2
)
(
2
1+X
,
ΨH (ζ ) = 2 ln
2
if it is unstable (ζ < 0). In the above,

X = (1 − 16ζ ) 1/4 .

(14.93)(14.94)(14.95)Using these values, convert the wind speed to that at 10 m in the neutral stability, and convert the temperature andspecific humidity to those at a height where the wind speed is observed,√[] ) −1
(
CDN10
zu
⃗
ln
− ΨM (ζu ),U10N = |∆U (zu )| 1 +
κ
z10
[
]
θ∗
zθ
θ(zu ) = θ(z θ ) −ln+ ΨH (ζu ) − ΨH (ζ θ ) ,
κ
zu
[
]
zq
q∗
q(zu ) = q(zq ) −ln+ ΨH (ζu ) − ΨH (ζ q ) ,κzu– 164 –

(14.96)(14.97)(14.98)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes

where z10 means z = 10 m. Estimate the bulk coeﬃcient at 10 m in the neutral stability using U10N , and then obtain thebulk coeﬃcient at a height (zu ) where wind speed is observed,
CDN10()] 2 ,zu
1−1/2− ΨM (ζu )1 + CDN10 CDN10 ln
κ
z10
(
) 1/2
CDu
CHN10
CDN10
(
)] ,
CHu = [
1zu
−1/2
1 + CHN10 CDN10 ln− ΨH (ζu )
κ
z10
(
) 1/2
CDu
CEN10
CDN10
(
)] .
CEu = [
zu
1
−1/2
ln
1 + CEN10 CDN10− ΨH (ζu )κz10
CDu = [

(14.99)(14.100)(14.101)

Repeat the procedures to calculate the bulk coeﬃcients using these bulk coeﬃcients with temperature and specific humidity
at z = zu , and recalculate virtual temperature (14.87) and the scales for friction velocity, temperature, and specific humidity(Equations of (14.88), (14.89), and (14.90)).

14.9.3 Appendix of Kondo (1975): BULKKONDO2In Appendix 1-3 of Kondo (1975), approximate formulae of the bulk transfer coeﬃcients are given in support of practicalapplications. The BULKKONDO2 option realizes this in MRI.COM. In this scheme, providing that atmospheric elements areobserved at various heights, latent and sensible heat fluxes are estimated at a height where air temperature and humidity
are observed. Momentum fluxes are also estimated at this height. Wind speed and bulk transfer coeﬃcients should becorrected to ones at the height where the air temperature and humidity are observed. In the program of bulk.F90, users
may specify BULKITER option in configure.in to do this iteration. The observed heights are specified in namelist
nml_srfflx_dbase (see Table 14.18).
⃗ (zu )| in
The following procedure is executed starting from the observed wind speed relative to the surface current |∆Uunits of [m/s] at a height of zu [m]:1. Estimate the neutral bulk transfer coeﬃcient at 10 m assuming that the observed wind speed is at 10 m in the neutral⃗ (zu )|). When the 10 m wind speed is already known, the neutral bulk transfer coeﬃcient isstability (U10 = |∆Ugiven as follows:p
CDN10 = 10−3 {ad + bd U10d },(14.102)where ad , bd , and pd are nondimensional numbers depending on the wind speed (Table 14.14).
2. Recalculate the roughness length using this neutral bulk transfer coeﬃcient at 10 m and correct the height for windspeed once again:{}−1/2
z0 = exp ln z10 − κ CDN10,(14.103)⃗ (zu )| ln(10/z0 )/ ln(zu /z0 ),U10 = |∆U

(14.104)

where z10 means z = 10 m.
3. Re-estimate the neutral bulk transfer coeﬃcient at 10 m using (14.102).
4. Operations 1 through 3 above are repeated twice or three times to determine CDN10 and U10 .
5. Convert the bulk transfer coeﬃcient to that in the neutral stability at the observed height of temperature and humidity(zh ):
κ2(14.105)
CDNh = []2−1/2
κCDN10− ln(z10 /zh )CHNh,ENh =

1/2
κCDNh
1/2
κ CDN10/CHN10,EN10 + ln(zh /z10 )– 165 –.

(14.106)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.9 Bulk transfer coeﬃcient
The neutral bulk transfer coeﬃcients at 10 m for evaporative flux and sensible heat flux (CEN10 and CHN10 respectively) are given as follows:
CEN10 = 10−3 {ae + beU10e + ce (U10 − 8) 2 },p

−3

CHN10 = 10

(14.107)p{ah + bh U10h + ch (U10 − 8) 2 },(14.108)

where ae,h , be,h , ce,h , and pe,h are nondimensional numbers depending on wind speed (Table 14.14).
6. Estimate the bulk transfer coeﬃcient in which stability is considered using the neutral bulk transfer coeﬃcientobtained above.
Table14.14 Non-dimensional parameters used in calculating the bulk transfer coeﬃcient based on Kondo (1975)

U10(m · s−1 )0.3-2.2
2.2-5
5-8
8-25
25-50

ad

ae

ah

bd

be

bh

ce

ch

pd

pe

ph

0
0.771
0.867
1.2
0

0
0.969
1.18
1.196
1.68

0
0.927
1.15
1.17
1.652

1.08
0.08580.0667
0.025
0.073

1.23
0.05210.01
0.008−0.0161.185
0.05460.01
0.0075−0.01700
0−0.0004
00
0
0−0.00045
0

−0.15
1
1
1
1

−0.16
1
1
1
1

−0.1571
1
1
1

The bulk transfer coeﬃcient depending on the stability is estimated in the following operations. First, the stability of
the atmospheric boundary layer is defined as follows:s = s0where
s0 =⃗ (zh )| is approximated byand |∆U|s0 |,|s0 | + 0.01

(14.109)

Ts − θ a⃗ (zh )| 2 [1.0 + log10 (10/zh )]2 )|∆U.

⃗ (zh )| = |∆U⃗ (zu )| ln(zh /z0 )/ ln(zu /z0 ).
|∆U(14.110)(14.111)

The atmospheric boundary layer is unstable if s > 0 (Ts − θ a > 0) and stable if s < 0 (Ts − θ a < 0). Finally, the bulk
transfer coeﬃcient depending on the stability is given as follows: 0
CD = 
 CDNh {0.1 + 0.03s√+ 0.9 exp(4.8s)}
 CDNh (1.0 + 0.47 s)

s < −3.3−3.3 ≤ s < 00 ≤ s,(14.112) 0
CE = 
 CENh {0.1 + 0.03s√+ 0.9 exp(4.8s)}
 CENh (1.0 + 0.63 s)

s < −3.3−3.3 ≤ s < 00 ≤ s,(14.113) 0CHNh {0.1 + 0.03s + 0.9 exp(4.8s)}
CH = 
 C (1.0 + 0.63√ s) HNhs < −3.3−3.3 ≤ s < 00 ≤ s.

(14.114)14.9.4 Numerical implementation
Because the iteration loop of the bulk formulae is not vectorized, the count of iteration is set twice and the repeating
operation is explicitly written in the source programs.

14.9.5 Usage• BULKKONDO2 for Kondo (1975)• BULKNCAR for Large and Yeager (2004; 2009)• BULKITER for computing bulk-transfer coeﬃcient using the iterative method.
– 166 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes

14.10 Specification of input data
The attributes of external forcing data should be specified by namelist nml_force_data. Table 14.15 explains variables
that are contained in this namelist block.
Table14.15 Namelist nml_force_data related to surface forcing data.

variable name
name
txyu
l2d

units

logicalfile_data
file_data_grid
imfrc
jmfrc
num_elminterval
num_data_max
ifstartlrepeatlinterpilinearluniform
loffsetoffsetlfactorfactorldoubleiverbose
ldefinedsec

ifstart(6)logicallogicallogicallogicallogicallogicallogicaldescription
name of data object
grid point on which data should be placed
data is 2D or not
file that contains reference vales for surface restoring
forcingdata file name of grid
grid size of data in x direction
grid size of data in y direction
the number of elements contained in fileregular time interval of data
the number of record contained in the file[ymdhms] of the first record of the input filerepeat use of climatological data or notinterpolate horizontally or notinterpolation methoddata is horizontally uniform or notdata is oﬀset or not
oﬀsetting factor
data is multiplied or not
multiplicative factor
input data is double precision or notstandard output of progress
data object with this name is defined or not

usage.true. for surface data, default = .false.)
needed if linterp = .true.

−1 : monthly−999 : steady forcing
default = .false.
default = .false.
1 : linear, 2 : spline
default = .false.
default = .false.
default = .false.
default = .false.
1 : extensive, 0 : minimum
default = .true.

Format of data files that should be prepared is shown in the following. Note that one data file should contain only one
element.

– 167 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.11 Surface data objects
Format of surface forcing data (trcref(_surf)_conf%file_data)

integer(4), parameter :: num_data_max = 365*4, nu = 99
! "num_data_max" is the number of data
integer(4) :: imfrc, jmfrc! data size
character(128) :: file_data, file_data_gridreal(4) :: precip(imfrc,jmfrc,num_data_max)
real(8) :: alonf(imfrc), alatf(jmfrc)
logical :: linterp! main data
open (unit=nu,file=file_data,access=direct,recl=4*imfrc*jmfrc)
do n = 1, num_data_max
write(unit=nu,rec=n) precip(:,:,n)end doclose(nu)! longitude/latitude of the main data

if (linterp) then ! If input data is horizontally interpolated in the model.
open (unit=nu,file=file_grid)
write(nu) alonf, alatf
close(nu)end if

14.11 Surface data objects
External data either read from file or received from AGCM are stored in the data set object type type_surf_dbase(Table 14.16). This object contains a sub object meteor_alt of type type_surf_meteor_alt (Table 14.17) that storesthe measurement height of marine meteorological variables. The name (ReadData for Reanalysis or FromCGCM for
AGCM) of the object and the measurement height of marine meteorological variables must be specified by namelist
nml_srfflx_dbase (Table 14.18).
Air-sea fluxes computed on the basis of external data are stored in the data set object type type_surf_product_ao(Table 14.19). This object contains sub objects givn_ao and calc_ao of the object type type_surf_bulk_ao (Table
14.20). The object calc_ao stores elements related to the computation of fluxes based on the bulk formula and the object
givn_ao stores the corresponding elements received by other components, that is, atmospheric model. Similarly, air-ice
fluxes are stored in type_surf_product_ai (Table 14.21), which contains sub objects givn_ai and calc_ai of theobject type type_surf_bulk_ai (Table 14.22).
Table14.16 Contents of the object type_surf_dbase related to external surface forcing data.

variable namenamemeteor_alt
x_windy_windtemp_airtdew_airsphm_airs_windunits
type_surf_meteor_alt
cm sec−1 or dyn cm−2
cm sec−1 or dyn cm−2
◦C◦C1
cm sec−1

description
name of data object
Observation altitude of this data set
X-ward wind or wind stressY-ward wind or wind stresssurface air temperature
dew point temperature
specific humidity
scalar surface windContinued on next page

– 168 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes
Table 14.16 – continued from previous page
units
descriptionhPa
sea level pressure
erg s−1 cm−2
upward or downward shortwave, summed
erg s−1 cm−2
upward or downward shortwave, visible, beam
erg s−1 cm−2
upward or downward shortwave, visible, diﬀuse
erg s−1 cm−2
upward or downward shortwave, near Infra-red, beam
erg s−1 cm−2
upward or downward shortwave, near Infra-red, diﬀuse
erg s−1 cm−2
upward or downward longwave◦C
sea surface temperature
cm s−1
precipitation water flux
cm s−1river discharge water flux
−1−2
erg s cmriver discharge heat flux
cm s−1snow fall water flux
cm s−1iceberg discharge water flux
erg s−1 cm−2
iceberg discharge heat flux
1
sea ice fraction

variable nameslpress
sw_flxsw_flx_visb
sw_flx_visd
sw_flx_nirb
sw_flx_nird
lw_flxsst
precipriver_w
river_hsnowfalliceberg_w
iceberg_h
frac_ice

Table14.17 Contents of the object type_surf_meteor_alt that stores measurement height of meteorological elements of external data.

variable namealt_wind
alt_temp
alt_sphm

units
m
m
m

description
altitude of wind
altitude of temperature
altitude of specific humidity

Table14.18 Namelist nml_srfflx_dbase.

variable name
name

units

description
name of the external dataset

l_marine_meteor
alt_wind_meter
alt_temp_meter
alt_sphm_meterlogicalmmmread marine meteorological data or not
observed height of windobserved height of temperatureobserved height of specific humidity

usage
ReadData (Reanalysis) or FromCGCM (AGCM)(required)required
required if l_marine_meteor = .true.
required if l_marine_meteor = .true.
required if l_marine_meteor = .true.

Table14.19 Contents of the object type_surf_product_ao that stores air-sea fluxes based on type_surf_dbase.

variable namenamel_calc_wind_stressl_calc_bulk_fluxes
sw_nt_ao
sw_nt_visb_ao
sw_nt_visd_ao
sw_nt_nirb_ao
sw_nt_nird_ao
lw_nt_aotsfc_ocn
wind_speed
givn_aocalc_aounits
logicallogicalerg s−1 cm−2
erg s−1 cm−2
erg s−1 cm−2
erg s−1 cm−2
erg s−1 cm−2
erg s−1 cm−2◦C
cm sec−1
type_surf_bulk_ao
type_surf_bulk_ao

description
name of data object
calculate wind stress for this data set
calculate fluxes based on a bulk formula for this data set
net shortwave, summed
net shortwave, visible, beam
net shortwave, visible, diﬀuse
net shortwave, near Infra-red, beam
net shortwave, near Infra-red, diﬀuse
net longwave, air-ocean
sea surface (skin) temperaturescalar wind speed
given from other model or forcing data filecalculated using bulk formulae– 169 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.12 Technical details
Table14.20 Contents of the object type_surf_bulk_ao that stores elements related to computation of air-sea fluxes
based on bulk formula.

variable namenamex_rel_wind
y_rel_wind
x_stress
y_stress
frcvl
rel_windlatentsensible
evapo

units
m sec−1
m sec−1dyn cm−2
dyn cm−2
m sec−1
m sec−1
erg s−1 cm−2
erg s−1 cm−2
cm sec−1

description
name of data object
surface wind relative to surface current
surface wind relative to surface current
X-ward surface stressY-ward surface stressfriction velocity
surface wind relative to surface current
latent heat flux, air-ocean
sensible heat flux, air-ocean
evaporative water flux, air-ocean

Table14.21 Contents of the object type_surf_product_ai that stores air-ice fluxes based on type_surf_dbase.

variable namenamel_calc_wind_stressl_calc_bulk_fluxes
sw_nt_aisw_nt_insw_nt_visb_ai
sw_nt_visd_ai
sw_nt_nirb_ai
sw_nt_nird_ai
lw_nt_aitsfc_ice
wind_speed
givn_aicalc_aiunits
logicallogicalerg s cm−2
erg s cm−2
erg s cm−2
erg s cm−2
erg s cm−2
erg s cm−2
erg s cm−2◦Cm sec−1
type_surf_bulk_ai
type_surf_bulk_ai

description
name of data object
calculate wind stress for this data set
calculate fluxes based on a bulk formula for this data set
net shortwave, summed, absorbed at the ice surface
net shortwave, summed, penetrate into the ice interior
net shortwave, visible, beam
net shortwave, visible, diﬀuse
net shortwave, near Infra-red, beam
net shortwave, near Infra-red, diﬀuse
net longwave, air-ice
ice surface temperaturescalar wind speed
given from other model or forcing data filecalculated using bulk formulae

Table14.22 Contents of the object type_surf_bulk_ai that stores elements related to computation of air-ice fluxes
based on bulk formula.

variable namenamex_rel_wind
y_rel_wind
x_stress
y_stress
frcvl
rel_windlatentsensible
sublimctntmpcsntmpunits
m sec−1
m sec−1dyn cm−2
dyn cm−2
m sec−1
m sec−1
W m−2
W m−2
−1
m sec of sea waterm sec−1
m sec−1

description
name of data object
surface wind relative to ice drift
surface wind relative to ice drift
X-ward surface stressY-ward surface stressfriction velocity
surface wind relative to ice drift
latent heat flux, air-ice
sensible heat flux, air-ice
sublimation water flux, air-ice
bulk transfer coeﬃcient × scalar windbulk transfer coeﬃcient × scalar wind14.12 Technical details
14.12.1 Work flow in the main program
Work flow in the main program is as follows. Note that the sea ice model is part of the surface flux module (surfflux).
In the surface flux module, air-sea fluxes are calculated first, then they are partly modified by the presence of sea ice.

– 170 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes

ogcm__run (ogcm.F90)|+ LOOP_CGCM_SCUP|+ cgcm_scup_ctl__cpl_loop_start|||+ cgcm_scup__get_a2o(Update surface fluxes ’FromCGCM’)|||+ cgcm_scup__get_atmos_a2o (Update surface atmospheric state ’FromCGCM’)|+ LOOP_TIME_STEP|||+ part_1|||+ surfflux_ctl__set_surf_state
|
|
|
|
|
+ force__everytime(Update external forcing ’ReadData’)
|
|
|
|
|
+ solarangle__main(Compute solar insolation angle)|||+ surfflux_ctl__main|||+ surfflux__air_ocean(Determine Air-Sea flux,||
either ’ReadData’ or ’FromCGCM’)|||+ surfflux__seaice(Flux modification due to the presence of Sea Ice)|||
+ surfflux__oceanflux(Determine fluxes at the sea surface)|||+ solarangle__shrtwv_flux (Determine veritcal distribution|of solar insolation flux)|+ end LOOP_TIME_STEP|+ cgcm_scup_ctl__cpl_loop_end|||+ cgcm_scup__put_o2a(Transfer surface oceanic state to coupler)|+ end LOOP_CGCM_SCUP

14.12.2 Coupling with an atmospheric modela. Standard caseIn the coupled mode, the ocean model receives all kinds of surface fluxes from the atmospheric model and the surface
oceanic state such as sea surface temperature, surface oceanic current, and sea ice state, are transferred to the atmosphericmodel. The coupling interval is usually longer than the time step of the ocean model and the fluxes are kept fixed in theocean model during the coupling cycle. See Section 17.9.2 for details of flux adjustment for conserving heat and fresh
water fluxes by correcting errors caused by interpolation and sea ice evolution during the coupling loop.
Work flow is summarized as follows:
1. Getting fluxes from the atmospheric model at the beginning of a coupling cycle in cgcm_scup__get_a2o(cgcm_scup.F90).
2. Getting fluxes for this step from get_flux_a2o (get_fluxes.F90).
3. Sea surface fluxes are determined after solving the sea ice model
4. Solve time evolution of the oceanic state
5. Send oceanic state to the atmospheric model– 171 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201714.12 Technical details

b. Recalculating and replacing fluxes in the ocean model

In coupled mode (SCUPCGCM option) or in child mode of nesting (SUB option) in which fluxes are taken from the parent
model (NGETFLUX option which is often used when the parent model is coupled with the atmosphere), the received fluxes
may be replaced by the ones calculated on its own by using the surface atmospheric states taken from the parent component,
its surface oceanic states, and its bulk formula. User may also specify the grid points where fluxes are replaced. This
option is mainly intended to be applied to coastal region of a child model of nesting in which the fluxes received from theparent component may not fully reflect the fine scale surface state due to complex topography. Without an appropriate
feedback mechanism at the sea surface, a child model may sometimes get unstable. This runtime option is controlled bynamelist nml_cgcm_recalc and nml_cgcm_recalc_replace, which is explained in Tables 14.23 and 14.24.
Table14.23 Namelist nml_cgcm_recalc.

variable namel_recalc_flux

units
logicaldescription
Recalculating of fluxes

usage
default = .false.

Table14.24 Namelist nml_cgcm_recalc_replace.

variable namel_recalc_replace_all
file_recalc_replace

units
logicalcharacterdescription
Replacing fluxes for all grid points
File that contains mask that specify
grid points where fluxes are replaced.

usage
default = .false.
real(8) 2d array (1/0) written with direct access14.12.3 Properties of moist air
In MRI.COM, a set of formulae given by Gill (1982) is used to compute properties of moist air. Physical constants of sea
water is listed on Table 2.1 of Chapter 2.
a. Saturation specific humidity

We consider computing saturation specific humidity in an environment with the sea surface temperature t [◦C] and thepressure p [hPa].
The relation between the vapor pressure e [hPa] and specific humidity q is give by
e/p = q/(ϵ + (1 − ϵ )q),(14.115)

where ϵ is the molecular weight ratio between water vapor and air:
ϵ = mw /ma = 18.016/28.966 = 0.62197.

(14.116)This is solved for the specific humidity asq = ϵ e/(p − (1 − ϵ )e).

(14.117)

The relative humidity γ is the ratio of the mixing ratio of the mass of vapor to the mass of dry air to the saturated one.
The mixing ratio (r) of the mass of vapor to the mass of dry air is given by
r=

q
.
1−q

(14.118)

Thus the relative humidity is computed using the specific humidity q and the saturation specific humidity qs asγ=

q(1 − qs ).
qs (1 − q)(14.119)

The saturation vapor pressure esw [hPa] of pure water vapor over a plane water surface is given by
log10 esw (t) = (0.7859 + 0.03477 t)/(1 + 0.00412 t),
where t is temperature between ±40 ◦C.

– 172 –

(14.120)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 14

Sea surface fluxes′ of water vapor at saturation is not exactly e but is given by
In air, the partial pressure esw
sw

The value of f w is given by′esw= f w esw .

(14.121)f w = 1 + 10−6 p(4.5 + 0.0006 t 2 ),(14.122)

where p is the pressure (in units [hPa]).
The saturation vapor pressure over a salt solution is less than over fresh water. For sea water, the reduction is about 2 %′ ).
(a factor of 0.98 should be applied to esw
The saturation vapor pressure esi of pure water vapor over ice is given by
log10 esi (t) = log10 esw (t) + 0.00422 t.

(14.123)

The saturation partial pressure esi′ in moist air is f i times esi . Values of f i are given correct to 1 part in 104 by (14.122).
Thusesi′ = f w esi .
(14.124)b. Specific heat of air

7 (
8q )
R 1−q+
27ϵ= 1004.6 × (1 + 0.8735q) J kg−1 K−1(14.125)

L v = 2.5008 × 106 − 2.3 × 103 t J kg−1 .

(14.126)

cpa =c. Latent heat

The latent heat of vaporization is given by

The latent heat of sublimation is given by
L s = 2.839 × 106 − 3.6(t + 35) 2 J kg−1 .

(14.127)

d. Air density ( ρ a )The equation of state of moist air of temperature T [K] under the pressure P [Pa] is given by
ρa =P≡ P/RTv,
RT (1 − q + q/ϵ )(14.128)

whereis called the virtual temperature andTv ≡ T (1 − q + q/ϵ ) = T (1 + 0.6078q)(14.129)

R = 287.04 J kg−1 K−1(14.130)is the gas constant of dry air.

– 173 –

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 15

Turbulence Closure Models
The surface boundary layer is made turbulent by wind injecting momentum, so vertical mixing may be induced even if thestratification is stable. The stratification near the surface may be made gravitationally unstable by the surface buoyancy
loss, so the vertical convection may be induced. Though these phenomena occur in small scales, they are important for thelarge scale oceanic processes. However, these are neither expressed by the fundamental equation of the general circulation
models∗ , nor resolved by large scale models. Therefore, general circulation models express the eﬀects of small scale
turbulence by making some closure assumptions called turbulence closure in which mixing is expressed by synoptic states(velocity and temperature (salinity) solution of the model).
There are mainly two approaches to the turbulence closure problem: statistical closure models and empirical approaches.
MRI.COM provides several options of statistical closure models, which are explained in this chapter. Note that an empiricalapproach of K-profile parameterization proposed by Large et al. (1994) is most popular in the ocean-climate modelingcommunity. However, the current version of MRI.COM does not take this approach. In a future version, an interface to theCommunity Vertical Mixing Project (CVMix; Griﬃes et al., 2015) will be implemented to test the empirical approaches.
Three statistical turbulence closure models are supported: So called Mellor-Yamada model (Mellor and Yamada 1982,
Mellor and Blumberg 2004) presented in Section 15.2, Noh and Kim (1999) model presented in Section 15.3, and a twoequation turbulence closure model using generic length-scale (GLS) equation by Umlauf and Burchard (2003) presented
in Section 15.4. Table 15.1 summarizes features of turbulence closure models.
Table15.1 Features and equation numbers of turbulence closure models explained in this chapter.

model option
prognosticated variables
surface boundary condition for TKE
formula for the length scalestability function

Mellor and Yamada
MELYAMq2 /2 (15.42)
specified (15.45)(15.49), (15.48)
Solution of (15.39)

Noh and Kim
NOHKIM
E(15.57)
flux (15.61)(15.75)(15.72), (15.73)

Umlauf and Burchard
GLS
E (15.83), Ψ (15.84)
flux (15.93)(15.78)(15.99), (15.100)

15.1 Statistical closure model
In the statistical closure models, the physical properties in the basic equations of motion for a Boussinesq fluid are separated
into averaged components and perturbed components (Reynolds decomposition), and then the equations are time averaged.
The expressions for averaged velocity U, averaged pressure P, and averaged potential temperature Θ are
∂Ui
= 0,
∂ xi
DU j
∂
1 ∂P
ρ
∂ [ 1 ( ∂Uk ∂U j )]+ ϵ jkl f k Ul =(−⟨uk u j ⟩) −
− gj
+
2ν
+
,
Dt
∂ xk
ρ0 ∂ x jρ0 ∂ x k
2 ∂xj
∂ xk
DΘ
∂
∂ ( ∂Θ )=(−⟨uk θ⟩) +
κ
,
Dt
∂ xk
∂ xk ∂ xk(15.1)(15.2)(15.3)

where D( )/Dt ≡ Uk ∂( )/∂ x k + ∂( )/∂t, g j is the gravity vector, f k is the Coriolis vector, ϵ i jk is the alternating tensor,
ν is viscosity, and κ is diﬀusivity. Averaged quantities (resolved by the general circulation model) are represented bycapital letters, and turbulent components (unresolved by the general circulation model) are represented by lower-case∗

It could be expressed by relaxing the hydrostatic approximation but this is for future work.

– 175 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201715.2 Mellor and Yamada model
letters. The statistical averages of the turbulent components are represented by ⟨ ⟩. The equation for salinity is similar tothat for temperature (15.3). If density is calculated from temperature and salinity, the whole expression would becomecomplicated. Therefore, density is assumed to be a function only of temperature here (we set ρ = βθ, where β is thecoeﬃcient of thermal expansion). Though these equations basically correspond to the MRI.COM governing equations(2.16)-(2.20), they now have contributions from subgrid-scale processes (⟨uk u j ⟩ and ⟨uk θ⟩), which is the focus of thischapter.
The equations (15.1)-(15.3) show that the statistical averages of the turbulent components (⟨uk u j ⟩ and ⟨uk θ⟩), which
are called as the second order moments, are necessary to calculate the evolution of large scale quantities. However, they
are unknown and their evolution must be solved explicitly if we do not use a parameterization to presume them on thebasis of the large scale quantities. Following Kantha and Clayson (2000), the equations for the second order moments are
expressed as follows:D⟨ui u j ⟩∂⟨ui u j ⟩ ]∂ [1 [ ∂⟨pui ⟩ ∂⟨pu j ⟩ ]+⟨uk ui u j ⟩ − ν+++ f k (ϵ jkl ⟨ul ui ⟩ + ϵ ikl ⟨ul u j ⟩)
Dt
∂ xk
∂ xk
ρ0 ∂ x j∂ xi[⟨ p ( ∂ui ∂u j )⟩⟨ ∂ui ∂u j ⟩
∂U j
∂Ui ]
β
= −⟨uk ui ⟩− ⟨uk u j ⟩−(g j ⟨ui θ⟩ + gi ⟨u j θ⟩) +
+
− 2ν
,
∂ xk
∂ xk
ρ0
ρ0 ∂ x j ∂ x i∂ xk ∂ xk
D⟨u j θ⟩⟨ ∂θ ⟩⟨ ∂u j ⟩]∂ [1 ⟨ ∂p ⟩+⟨uk u j θ⟩ − κ u j−ν θ+θ
+ ϵ jkl f k ⟨ul θ⟩
Dt
∂ xk
∂ xk
∂ xk
ρ0 ∂ x j⟨ ∂u j ∂θ ⟩[∂U j ]β∂Θ− g j ⟨θ 2 ⟩ − (κ + ν),
− ⟨uk θ⟩= −⟨u j uk ⟩
∂ xk
∂ xk
ρ0
∂ xk ∂ xk⟨ ∂θ ∂θ ⟩⟨ ∂θ 2 ⟩]D⟨θ 2 ⟩∂ [∂Θ+
⟨uk θ 2 ⟩ − κ= −2⟨uk θ⟩
− 2κ
.
Dt
∂ xk
∂ xk
∂ xk
∂ xk ∂ xk(15.4)(15.5)(15.6)

Thus, to solve the evolution of the second order moments, information about higher moments such as ⟨uuu⟩, ⟨uuθ⟩, ⟨uθ 2 ⟩is necessary. The system of equations is not closed, and assumption must be made at some level to close the problem.
The expression closed in the second order of the turbulent components is named second moment closure and is frequently
used in modeling turbulence. Many models have been proposed to parameterize the higher order correlation terms in (15.4)
through (15.6) in order to close the problem at the second moment. Those models diﬀer mainly in the manner how thepressure-strain correlators are parameterized (Burchard and Bolding, 2001). After applying independent parameterizations
for the higher order terms, the majority of models adopt several similar simplifications such as algebraization and boundary
layer approximation that result in the "stability functions" of their own. Stability functions are determined by local shear
and stability and relate the vertical gradient of large scale fields with vertical turbulent fluxes. Procedures adopted bymajor statistical closure models are summarized and compared by Burchard and Bolding (2001).

15.2 Mellor and Yamada model
Mellor and Yamada (1982) model is a classic model as well as a milestone in the eﬀorts of finding practical solutions formarine and atmospheric boundary layers. We show some details of the simplifying procedures made in this model as a
typical example. This is also intended to be an introduction to the problem of turbulence closure. However, it should be
noted that this is neither unique nor final solution.

15.2.1 Fundamental closure assumptions
Mellor and Yamada (1982) closed the system (15.1)-(15.6) by reducing the higher order terms as follows. Based on Rotta’s(1951a,1951b) hypothesis of energy redistribution, the covariances of pressure and velocity gradients are assumed to belinear functions of Reynolds stress:⟨ p ( ∂ui
ρ0 ∂ x j+( ∂Ui ∂U j )∂u j )⟩
δi j 2 )q (=−⟨ui u j ⟩ −q + C1 q2
+
,
∂ xi
3l 1
3
∂xj
∂ xi

(15.7)

where q2 ≡ ⟨ui2 ⟩, l 1 is the length scale, C1 is a non-dimensional constant, and δi j is Kronecker’s delta, which is unity fori = j and zero for i , j.
Using Kolmogorov’s hypothesis of local isotropy in small eddies, the energy dissipation is modeled as follows:2ν⟨ ∂ui ∂u j ⟩∂ xk ∂ xk=−

where Λ1 is the length scale.
– 176 –2 q3
δi j ,
3 Λ1(15.8)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 15

Turbulence Closure Models

The redistribution of temperature and the dissipation of heat are modeled in the same form as above:⟨ p ∂θ ⟩q= − ⟨u j θ⟩
ρ0 ∂ x j
3l 2(κ + ν)⟨ ∂u j ∂θ ⟩∂ xk ∂ xk= 0,(15.9)(15.10)where l 2 is the length scale. The dissipation of temperature variance is2κ⟨ ∂θ ∂θ ⟩q= −2 ⟨θ 2 ⟩,∂ xk ∂ xk
Λ2(15.11)

where Λ2 is the length scale.
In order to avoid the higher order problems, the turbulent velocity diﬀusion term and the other higher order terms aremodeled as follows:( ∂⟨ui u j ⟩ ∂⟨ui uk ⟩ ∂⟨u j uk ⟩ )3++,⟨uk ui u j ⟩ = lqSq
5
∂ xk
∂xj
∂ xi
( ∂⟨uk θ⟩ ∂⟨u j θ⟩ )⟨uk u j θ⟩ = − lqSuθ
,
+
∂xj
∂ xk
⟨uk θ 2 ⟩ = − lqSθ∂⟨θ 2 ⟩,
∂ xk

(15.12)(15.13)(15.14)

where Sq , Suθ , and Sθ are non-dimensional numbers and can be set as constants or functions of certain parameters. Other
relations are ⟨pθ⟩ = 0 and ⟨pui ⟩ = 0.
The essence of the Mellor-Yamada mixed layer model is that the above length scales are related linearly to each other:(l 1, Λ1, l 2, Λ2 ) = ( A1, B1, A2, B2 )l,(15.15)where l is the vertical scale of turbulence (also called the master length scale), and A1, B1, A2, B2 , and C1 are empirical constants and are determined from experiment data. Mellor and Yamada (1982) employ ( A1, B1, A2, B2, C1 ) =(0.92, 16.6, 0.74, 10.1, 0.08).

15.2.2 Algebraization, boundary layer approximation, and stability functions for the level 2.5 model
The turbulence model that solves the evolution of the statistically averaged values of the second-order turbulent components
based on the simplification described in the previous subsection is called the level-4 model.
The level-3 model solves the evolution of the turbulent kinetic energy (q2 /2) and the variance of potential temperature(⟨θ 2 ⟩) (in some cases, the covariance of potential temperature and salinity (⟨θs⟩) and the variance of salinity (⟨s2 ⟩)). The
other statistically averaged values are solved diagnostically through algebraic equations assuming them to be in the steadystate.
In the level-2.5 model, the variance of the potential temperature is also assumed to be in a statistically steady state (seeexpression (15.33) that appears later).
In the level-2 model, the turbulent kinetic energy is also assumed to be in a statistically steady state.
The level-2.5 model, which is employed as a surface boundary layer model by MRI.COM, is further simplified byapplying the following boundary layer approximations.
• Neglect the Coriolis term in the equations of motion for the turbulent components.
• Neglect the molecular viscosity and diﬀusivity.
• Use the hydrostatic assumption in the vertical component of the equation of motion.
• Consider only vertical diﬀerentiation (direction perpendicular to the boundary) in the spatial diﬀerentiation for theterm involving turbulent velocity.

– 177 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201715.2 Mellor and Yamada model
The equations for the large scale quantities, (15.2) and (15.3), become
DU
∂
1 ∂P
+
⟨uw⟩ = −+ f V,
Dt∂zρ0 ∂ xDV∂
1 ∂P+⟨vw⟩ = −− f U,
Dt∂zρ0 ∂ y
ρ
1 ∂P
−g ,
0= −
ρ0 ∂z
ρ0
DΘ
∂
1 ∂R
+
(⟨wθ⟩) =.
Dt∂zρ0 cp ∂z(15.16)(15.17)(15.18)(15.19)In (15.19), R is the heat flux due to penetration of shortwave radiation. The level-2.5 system consists of the time evolution
equation for turbulent kinetic energy and algebraic equations for other second-moment turbulent quantities.
The time evolution equation for the turbulent kinetic energy is
D ( q2 )∂ [∂ ( q2 )]−lqSq= Ps + Pb − ϵ,
Dt 2
∂z
∂z 2
where
Ps = −⟨wu⟩∂U∂V− ⟨wv⟩∂z∂z(15.20)(15.21)is the term for energy produced by the vertical shear of the averaged flow,
Pb = −g⟨w ρ⟩/ρ0(15.22)ϵ = q3 /Λ1(15.23)is the term for energy produced by buoyancy, andis the energy dissipation term. The constant Sq is given in the next subsection.
The algebraic equations for the statistically averaged values, which are expressed by other second-moment turbulent
quantities, are given below.
]
q2 l 1 [∂U∂V+−4⟨wu⟩+ 2⟨wv⟩− 2Pb ,
3
q
∂z
∂z
2
[
]
q
l1
∂U
∂V
⟨v 2 ⟩ =+2⟨wu⟩− 4⟨wv⟩− 2Pb ,
3
q
∂z
∂z
2
[
]
q
l1
∂U
∂V
⟨w 2 ⟩ =+2⟨wu⟩+ 2⟨wv⟩+ 4Pb ,
3
q
∂z
∂z
∂V
∂U ]
3l 1 [⟨uv⟩ =−⟨uw⟩− ⟨vw⟩,q∂z∂z]
3l 1 [∂U⟨wu⟩ =−(⟨w 2 ⟩ − C1 q2 )− g⟨u ρ⟩ ,
q
∂z
[
]
∂V
3l 1
⟨vw⟩ =−(⟨w 2 ⟩ − C1 q2 )− g⟨v ρ⟩ ,q∂z
3l 2 [∂Θ∂U ]⟨uθ⟩ =−⟨uw⟩− ⟨wθ⟩,q∂z∂z
3l 2 [∂Θ∂V ]⟨vθ⟩ =−⟨vw⟩− ⟨wθ⟩,q∂z∂z]
3l 2 [∂Θ⟨wθ⟩ =−⟨w 2 ⟩− g⟨θ ρ⟩ ,
q
∂z
Λ
∂Θ
2
⟨θ 2 ⟩ = − ⟨wθ⟩.
q∂z
⟨u2 ⟩ =

(15.24)(15.25)(15.26)(15.27)(15.28)(15.29)(15.30)(15.31)(15.32)(15.33)

Some of the terms in these equations, which are used in (15.16)-(15.19), can be further deformed by substitution as– 178 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 15

Turbulence Closure Models

follows:∂U,∂z∂V−⟨vw⟩ = K M,∂z∂Θ−⟨θw⟩ = K H,∂z
K M = lqSM ,K H = lqSH .
−⟨uw⟩ = K M(15.34)(15.35)(15.36)(15.37)(15.38)This simplification means that the vertical turbulent fluxes are proportional to the gradient of the large scale fields. The
ultimate purpose of solving the mixed layer model is to determine the coeﬃcients of momentum and heat fluxes, K M andK H , using (15.37) and (15.38).
Assuming that the potential density is linearly related to the potential temperature (and salinity), the simultaneous
equations for SM and SH are derived as follows:
SM [6A1 A2 g M ] + SH [1 − 3A2 B2 gH − 12A1 A2 gH ] = A2,
SM [1 + 6A21 g M − 9A1 A2 gH ] − SH [12A21 gH + 9A1 A2 gH ] = A1 (1 − 3C1 ),(15.39)wherel 2 [( ∂U ) 2 ( ∂V ) 2 ]l2+= 2 S 2,
2
∂z
∂z
q
q
2
2
l g ∂ ρ̃lgH ≡ 2= 2 N 2,q ρ0 ∂zq

gM ≡(15.40)(15.41)and ∂ ρ̃/∂z is the vertical gradient of potential density. Note that SM and SH are stability functions of this closure model.
Using SM and SH , K M and K H are then obtained from (15.37) and (15.38) by determining q and l.

15.2.3 Solving for velocity scale (q) and length scale (l )The turbulent velocity q is obtained by solving the following expression that is modified from (15.20) using the above
results:[ ( ∂U ) 2 ( ∂V ) 2 ]∂ ( q2 )∂ [g∂ ( q2 )]
∂ ρ̃
−
KE
= KM
+
+
KH
− ϵ,
(15.42)
∂t 2
∂z
∂z 2
∂z
∂z
ρ0
∂z
where
K E = lqSq

(15.43)and advection terms are neglected. In MRI.COM, Sq is set proportional to SM (Sq ∝ SM ). We adopt the form
Sq = Sqc SM /SM n , where Sqc = 0.2 and SM n = 0.3927. With this choice, Sq = 0.2 when the stratification is neutral(gH = 0).
The sea surface boundary condition for the turbulent kinetic energy follows Mellor and Blumberg (2004):
Kq∂q2= 2αCB uτ3,∂z(15.44)

where αCB = 100 and uτ is the frictional velocity defined as uτ ≡ (τs /ρs ) 1/2 by using the surface stress (τs ) and the sea
surface density (ρs ). Mellor and Blumberg (2004) showed that this flux boundary condition is analytically converted tothe condition for q at the sea surface:q2 = (15.8αCB ) 2/3 uτ2 .
(15.45)On the other hand, the bottom boundary condition isq = 0.

(15.46)The vertical scale of the turbulence (master length scale, l) is estimated by many formulae such as a time evolution
equation (which is usually empirical and is not completely based on physics) and a diagnosis. MRI.COM uses diﬀerent
diagnosis equations for the surface layer influenced directly by surface wind wave eﬀects and for the internal region. The– 179 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201715.2 Mellor and Yamada modelsurface layer is defined by a depth of |z| < zw , where the roughness parameter zw due to surface wind waves is given by
Mellor and Blumberg (2004) as follows:u2z w = βw τ ,gβw = 2.0 × 105 .

(15.47)In this layer, the formula for l is
l = κzw,(15.48)where κ is the von Karman constant (κ = 0.4). In the internal region (|z| > zw ), the formula for l is
l=γ

∫ 0

′

|z |qdz

′

/∫ 0

zb

qdz ′,(15.49)zb

where γ = 0.2, and z b is the water depth. This is recognized as the averaged depth with the weight of the kinetic energy,
which is suﬃcient for the ocean boundary layer according to Mellor and Yamada (1982).

15.2.4 Implementation
This section briefly describes the solution procedure.
The mixed layer model (subroutine name mysl25 in my25.F90) is called as the last procedure of each time step thatproceeds from n to n+1. After the master length scale (l) for the present time step (n) is determined using (15.48) and(15.49), where the latter is actually computed in the previous time step, the turbulent kinetic energy (q2 /2) is solved using(15.42) with (15.23) under the boundary conditions of (15.45) and (15.46), where the forward finite diﬀerence (n → n+1)
is used in time. The implicit method is used for the vertical diﬀusion of the turbulent kinetic energy and energy dissipation
term, since these terms could become significantly large (see Section 19.5). The vertical viscosity and diﬀusivity for thetime step n+1 are estimated using q, l, and (15.37) to (15.39). These are used in the governing equations for the synoptic
scale field. The vertical scale of the turbulence (master length scale) based on the turbulence kinetic energy (15.49) is
calculated to prepare for the next time step.
The turbulent kinetic energy and the master length scale are defined at the bottom of the tracer cell (i, j, k). The specificexpression for the discretized form of the turbulent kinetic energy (E = q2 /2) equation, (15.42), is as follows:
Ekn+1 − Ekn
∆t

=

n+1
n+1
n+1
n+1
1 [ K E k− 12 (Ek−1 − Ek ) K E k+ 21 (Ek − Ek+1 ) ]−∆zk∆zk− 1∆zk+ 1
2+ KM k− KH k2(U n+11 − U n+11 )(Ũk− 1 − Ũk+ 1 )
k− 2

k+ 2

2

2

∆zk2
B n+11 − B n+11
k− 2

∆zk

k+ 2

+ KM k

(V n+11 − V n+11 )(Ṽk− 1 − Ṽk+ 1 )
k− 2

k+ 2

∆zk2

2

2

(15.50)− 2Ekn+1 qkn /B1 l kn,

where Ũ = (U n+1 + U n )/2, Ṽ = (V n+1 + V n )/2, and B is buoyancy (= − gρρ0 ). The discrete expression for shear production(the second and third terms on the r.h.s.) and the buoyancy sink (the fourth term on the r.h.s.) follows Burchard (2002),
which is consistent with the conservation law of the sum of mean and turbulent energy.
To summarize, the numerical operations proceed in the following order:
1. Calculate the master length scale for the present time step using (15.47) and (15.48) as well as (15.49) from theprevious time step.
2. Update the turbulent kinetic energy using (15.50), (15.45) and (15.46).
3. Solve the algebraic equation for SM and SH using (15.39).
4. Calculate the vertical viscosity and diﬀusivity for the next time step using (15.37), (15.38), and (15.43).
5. Calculate the master length scale using (15.49) in preparation for the next time step.

15.2.5 Usage
Mellor and Yamada model is invoked by MELYAM option on compilation. There is no tuning parameter at run time. Initialstate for vertical diﬀusion of turbulence kinetic energy (K E ), turbulence kinetic energy (q), and length scale (l) is readfrom restart file. Otherwise it is set as a state without turbulence. The state without turbulence is zero turbulence kinetic– 180 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 15

Turbulence Closure Models

energy (q = 0), the mixing length of 1 meter (l = 1 m) and the vertical diﬀusion of turbulence kinetic energy is the constantbackground vertical viscosity given by visc_vert_bg_cm2ps in namelist nml_visc_vert_bg. How model is initializedis specified by namelist nml_melyam_run. Parameters are listed on Table 15.2.
Table15.2 namelist nml_melyam_run

variable namel_rst_melyam_in

units
logicaldescription.true. : Read restart files specified by nmlrs_my_avq,
nmlrs_my_q, and nmlrs_my_alo for the initial condition.
.false.: Start condition is that of no turbulence. See text
for details.

usage
default = l_rst_in of
nml_run_ini_state

15.3 Noh and Kim (1999) model
The mixed layer model proposed by Mellor and Yamada was originally developed for the atmospheric boundary layer, andits surface boundary is treated as a solid wall. When they applied this model to the ocean, they regarded that the turbulent
kinetic energy is injected into the ocean by the wind stress at the solid-wall sea surface. The model by Mellor and Yamada
could therefore be considered to insuﬃciently represent the oceanic turbulent mixed layer.
Noh and Kim (1999) presented a model that can resolve this insuﬃciency. The model is basically the same as Mellor
and Yamada level 2.5 scheme in that turbulent kinetic energy is prognosticated and the length scale is diagnosticallydetermined. It diﬀers in the stability function and the treatment of surface boundary condition.

15.3.1 Fundamental equation
The equations for the zonal and meridional components of the velocity, U, V , buoyancy B = −g∆ρ/ρ0 , and turbulent
energy E in the large scale fields are
DU
∂
1 ∂P
= − ⟨uw⟩ −+ f V,
Dt∂zρ0 ∂ xDV∂
1 ∂P= − ⟨vw⟩ −− f U,
Dt∂zρ0 ∂ yDB∂∂R= − ⟨bw⟩ +
,
Dt
∂z
∂z
)⟩
DE
∂U
∂ ⟨ ( p∂V=−w
+ uu + vv + ww − ⟨uw⟩
− ⟨vw⟩+ ⟨bw⟩ − ϵ,
Dt
∂z
ρ0
∂z
∂z

(15.51)(15.52)(15.53)(15.54)

where R is the downward shortwave radiation and ∂ R/∂z is its convergence.
The turbulent flux is expressed by using the large scale fields (in capital letters) as follows:
DU
∂ (
∂U )
1 ∂P
=
KM
−
+ f V,Dt∂z∂zρ0 ∂ x
DV
∂ (
∂V )
1 ∂P
=
KM
−
− f U,Dt∂z∂zρ0 ∂ yDB∂ (∂B ) ∂ R
=
KH
−
,
Dt
∂z
∂z
∂z
DE
∂ (
∂E )
∂U ∂U
∂V ∂V (
∂B )
=
KE
+ KM
+ KM
− KH
− ϵ.
Dt
∂z
∂z
∂z ∂z
∂z ∂z
∂z

(15.55)(15.56)(15.57)(15.58)The surface boundary conditions are as follows:
∂U
τ
= ,
∂z
ρ0
∂B
KH
= Q0,
∂z
∂E
KE
= muτ3,∂z

KM– 181 –

(15.59)(15.60)(15.61)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201715.3 Noh and Kim (1999) model
where m is a tuning parameter and m = 100 is recommended by Noh and Kim (1999). Unlike Mellor and Yamada model,
a flux boundary condition is used for turbulence kinetic energy. This is intended to express the direct input of turbulencekinetic energy due to wind waves. As such, in the presence of sea ice, the direct input of wave energy may be more
appropriate only through open water. The boundary condition may be optionally modified so thatKE∂E
= muτ3 (1 − a),∂z(15.62)where a is area fraction of sea ice. The no-flux condition is used for the bottom:
KE

∂E
= 0.
∂z

(15.63)15.3.2 Stability function
The central problem is how to determine the viscosity, diﬀusivity (K M , K H , K E ), and turbulent energy dissipation rate (ϵ).
In the Noh and Kim (1999) scheme, by using the typical velocity scale (q = (2E) 1/2 ) and the vertical length scale (l) of
the turbulence, they are obtained by the following formula:
K M = Sql,
K H = SB ql,
K E = SE ql,
ϵ = Cq3 l −1 .

(15.64)(15.65)(15.66)(15.67)The constants (S, SB, SE , C) are obtained from experiments. For neutral stratification, it is assumed that S = 0.39,
Pr ≡ S/SB = 0.8, σ ≡ S/SE = 1.95, and C = 0.06. Hereafter, these values for neutral stratification are denoted as S0 ,Pr0 , σ0 , and C0 , respectively.
Regarding the influence of the stratification, that is, stability function, we assume that the vertical scale of turbulence islimited by the vertical scale of buoyancy l b = q/N (N 2 = ∂B/∂z). That is,where Rit is the turbulent Richardson numberK ∼ ql b ∼ qlRit −1/2,(15.68)

Rit = (Nl/q) 2 .

(15.69)Rit corresponds to gH of Mellor and Yamada model. See (15.41).
This means that when the stratification is strong (N is large, Rit is large, and K is small), the turbulent energy is nottransported downwards. It could also be considered that the local turbulent energy dissipation becomes large.
The following equation is used for S so that it satisfies (15.68) when Rit is large:
S/S0 = (1 + αRit ) −1/2,(15.70)where α is a tuning parameter. Noh and Kim (1999) recommend α ∼ 120.0, but α ∼ 5.0 is the default value of MRI.COM.
The eﬀect of stratification on the Prandtl number (Pr) is set following Noh et al. (2005):
Pr/Pr0 = (1 + βRit ) 1/2,

(15.71)

where β is a tuning parameter and 0.5 is used following Noh et al. (2005).
To summarize, the stability function for this model is expressed as follows:
S0(1 + αRit ) 1/2
S
S0
SB =
=
.
Pr Pr0 (1 + αRit ) 1/2 (1 + βRit ) 1/2
S=(15.72)(15.73)The eﬀect of stratification on the energy dissipation (C) is set as follows:
C/C0 = (1 + αRit ) 1/2 .
Note that this is independent of stratification for Mellor and Yamada and Umlauf and Burchard (2003) model.
– 182 –

(15.74)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 15

Turbulence Closure Models

In the case of unstable stratification (N 2 < 0), K M = K H = 1.0 m2 s−1 and K E is estimated from the turbulent velocity
scale and the vertical length scale in the model. This treatment is due to the diﬀerence between the time scales of thevertical convection and the development of turbulence.
Finally, the vertical scale of turbulence is diagnosed asl=κ(|z| + z0 ),(1 + κ(|z| + z0 )/h)(15.75)where z0 is the sea surface roughness (z0 = 1 m), z is the depth, and h is the mixed layer depth. The mixed layer depth (h)is determined as the depth at which the buoyancy frequency becomes the maximum in the vertical column. The verticalscale becomes longer as the mixed layer becomes deeper.

15.3.3 Implementation
Equation (15.58) is solved for the prognostic variable E in the subroutine nkoblm in nkoblm.F90 as the last procedure ofeach time step. The forward finite diﬀerence is used in the time evolution. The implicit method is used for the vertical
diﬀusion of the turbulent kinetic energy and energy dissipation term, since these terms could become significantly large(see Section 19.5). The new E is used to determine the coeﬃcients of viscosity and diﬀusivity for the next time step.
The turbulent kinetic energy and the master length scale are defined at the center of the tracer cell (i, j, k − 12 ). The
specific expression for the discretized form of the turbulent kinetic energy (E) equation is as follows:
E n+11 − E n 1
k− 2

k− 2

∆t

1 [
=
∆zk− 1K E k−1 (E n+13 − E n+11 )2

1+ K M k−1
2
1
+ K M k−1
2
1
− K H k−1
2

k− 2

k− 2

∆zk−1

−

K E k (E n+11 − E n+11 ) ]
k− 2

∆zk

(un+13 − un+11 )(un+13 − un 1 )
k− 2
k− 2
k− 2
k− 2
∆zk−1 ∆zk− 1
2(v n+13 − v n+11 )(v n+13 − v n 1 )
k− 2
k− 2
k− 2
k− 2
∆zk−1 ∆zk− 1
2

B n+13 − B n+11
k− 2

∆zk−1

k− 2

1
− KH k2

k+ 2

1+ KM k2
1+ KM k2

B n+11 − B n+11
k− 2

∆zk

k+ 2

(un+11 − un+11 )(un 1 − un+11 )
k− 2

k+ 2

k− 2

k+ 2

∆zk ∆zk− 1
2(v n+11 − v n+11 )(v n 1 − v n+11 )
k+ 2
k− 2
k+ 2
k− 2
∆zk ∆zk− 1(15.76)2n+1 n/l n .
− 2CEk−1q
k− 1 k− 12

2

2

The discrete expression for the shear production (the second through fifth terms on the r.h.s.) and the buoyancy sink(the sixth term on the r.h.s.) follows Burchard (2002), which is consistent with the conservation law of the sum of mean
and turbulent energy.
To summarize, the numerical operations proceed in the following order:
1. Update the master length scale using (15.75),
2. Update the turbulent kinetic energy using (15.76) under the boundary conditions (15.61) and (15.63)
3. Calculate the vertical viscosity and diﬀusivity for the next time step (15.64) and (15.65).

15.3.4 Usage
Noh and Kim (1999) model is invoked by NOHKIM option on compilation.
Initial state for vertical diﬀusion of turbulence kinetic energy (K E ), turbulence kinetic energy (E) is either read fromrestart file or set as a state without turbulence. The state without turbulence means that E is very small as 12 qmin 2 , whereqmin = 1 × 10−5 m s−1 is the minimum of turbulence velocity scale, and that K E is set the constant background verticalviscosity given by visc_vert_bg_cm2ps in namelist nml_visc_vert_bg. How model is initialized is specified bynamelist nml_nohkim_run. Parameters are listed on Table 15.4.
Tuning parameters may be specified at run time by namelist nml_nkoblm. Parameters are listed on Table 15.3.

– 183 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201715.4 Generic length scale model by Umlauf and Burchard (2003)
Table15.3 namelist nml_nkoblm

variable name
mfnk
alpnk
betankl_wave_damp_seaice

units
1
1
1
logicaldescription
m of Eq. (15.61)
α of Eq. (15.70)
β of Eq. (15.71)
Eq. (15.62) is used instead of (15.61)

usagedefault = 100.0
default = 5.0
default = 0.5
default = .false.

Table15.4 namelist nml_nohkim_run

variable namel_rst_nohkim_in

units
logicaldescription.true. : Read restart files specified by nmlrs_nk_avq
and nmlrs_nk_eb for the initial condition.
.false.: Start condition is that of no turbulence. See text
for details.

usage
default = l_rst_in of
nml_run_ini_state

15.4 Generic length scale model by Umlauf and Burchard (2003)†

The turbulent closure models of Mellor and Yamada level 2.5 (section 15.2) and Noh and Kim (section 15.3) asincorporated in MRI.COM have only one prognostic variable of turbulence (turbulent kinetic energy). These may beclassed as "one-equation" turbulent closure models. In these models, some formulae for length scale are given (e.g.,(15.49) and (15.75)). Therefore the turbulent length scale l is not a predicted but a diagnosed quantity.
Generic length scale model by Umlauf and Burchard (2003) is a generalized form of two-equation turbulence models.
It introduces a generic length-scale variable ψ as the second variable, expressed as
ψ = (cµ0 ) p E m l n .

(15.77)

With E and ψ given, the turbulent length scale can be computed aspl = (cµ0 ) − n E − n ψ n ,m

1(15.78)and the dissipation rate ϵ can be computed from the following relation,
ϵ = (cµ0 ) 3

E 3/2
.
l

(15.79)With appropriate values of the exponents p, m, n in (15.77), a number of well-known two-equation models (e.g., k-ϵ modeland k-ω model) can be identified as special case of the generic model.

15.4.1 Fundamental equation
The fundamental equations for U, V , B and E are basically common to those of Noh and Kim model. The equation forgeneric length scale variable ψ is added. Expressing the turbulent flux using the large scale fields with the assumption ofhorizontally homogeneous flow, the fundamental equations of the model is expressed as follows:
DU
∂ (
∂U )
1 ∂P
=
KM
−
+ f V,Dt∂z∂zρ0 ∂ x
DV
∂ (
∂V )
1 ∂P
=
KM
−
− f U,Dt∂z∂zρ0 ∂ y∂ (∂B ) ∂ R
DB
=
KH
−
,
Dt
∂z
∂z
∂z
DE
∂E )
∂ (
=
KE
+ Ps + Pb − ϵ,
Dt
∂z
∂z
Dψ
∂ ( ∂ψ ) ψ=Kψ+ (cψ1 Ps + cψ3 Pb − cψ2 ϵ ),
Dt
∂z
∂z
E
† This model is implemented by Shiro Nishikawa as he was a visiting researcher at MRI.

– 184 –

(15.80)(15.81)(15.82)(15.83)(15.84)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 15

Turbulence Closure Models

where Ps and Pb are the shear production and the buoyancy production expressed according to
Ps = −⟨uw⟩∂u∂v− ⟨vw⟩= K M M 2,
∂z
∂z

and

Pb = ⟨bw⟩ = −K H N 2,

(15.85)( ∂v ) 2∂B,
N2 =(15.86)∂z∂z∂z
are the shear frequency and the buoyancy frequency, respectively. The vertical turbulent diﬀusivities in horizontallyhomogeneous flows are expressed asM2 =( ∂u ) 2+

K M = cµ E 1/2 l,(15.87)

K H = cµ′ E 1/2 l,(15.88)ψl/σk ,(15.89)Kψ = cµ E 1/2 l/σψ,(15.90)K E = cµ E

1/2ψwhere cµ and cµ′ are stability functions and σk and σψ are the constant Schmidt numbers. Here, cµ and cµ′ corresponds to
SM , and SH of the Mellor and Yamada model with a relation√√cµ = 2SM , cµ′ = 2SH .
(15.91)In the case of unstable stratification (N 2 < 0), K M = K H = 1.0 m2 s−1 is used in MRI.COM. This treatment is common
to Noh and Kim model (section 15.3).

15.4.2 Boundary conditions
At the surface, the shear-free boundary layers with injection of turbulent kinetic energy are assumed, in which
E = K (−z + z0 ) α,l = L(−z + z0 ),(15.92)are supposed, where K is the scale of turbulent kinetic energy, L and α are empirical constants given in Table 15.5, andz0 is the sea surface roughness (z0 = 1 m). K is determined by the injection of turbulent kinetic energy which is given by
the surface flux. We use flux forms for the surface boundary conditions (z = 0) of E and ψ as follows:∂E= ηuτ3,∂zcµ (cµ0 ) p
1
1
∂ψ
=−
(mα + n)K m+ 2 L n+1 (z0 ) (m+ 2 )α+n .
Fψ = Kψ∂zσψFE = K E(15.93)(15.94)

As explained in the previous section for Noh and Kim (1999) model, injection of turbulence kinetic energy may be madeonly through the open water in the presence of sea ice.
FE = K E∂E= ηuτ3 (1 − a),∂z(15.95)where a is sea ice area fraction.
The parameter K can be determined from
ψ

σk

2
3

1
K = *−
FE + α ,
c
αL
z
, µ
- 0

(15.96)

where (15.89), (15.92) and (15.93) are used.
At the bottom, the logarithmic boundary layers are assumed, in which l = κ(z + H + z0b ) is supposed. The bottom
boundary conditions (z = −H) for E and ψ are given by flux forms as follows:
∂E
= 0,
∂z
0 p+1 n+11
∂ψ n(cµ ) κFψ = Kψ=E m+ 2 (z0b ) n,∂zσψFE = K E(15.97)(15.98)where z0b is the bottom roughness, which is set to be z0b = 1 m. Implementation of the bottom boundary condition isrealistic as compared to the Mellor and Yamada scheme and the Noh and Kim scheme.
– 185 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201715.4 Generic length scale model by Umlauf and Burchard (2003)

15.4.3 Model parameters
Umlauf and Burchard (2003) investigated the properties of the generic model in some fundamental flows and derived theψconstraints on the model parameters (cµ0 , cψ1 , cψ2 , cψ3 , σk , σψ , m and n). Table 15.5 is a recommended parameter set
for the generic model by Umlauf and Burchard (2003), based on their calibration. MRI.COM uses this as the default
parameters. Note that the parameter p (the factor (cµ0 ) p in (15.77)) is mathematically irrelevant, since (15.84) can be
multiplied by any constant without changing the solution (p = 2.0 is used as the default value).
Table15.5 A recommended parameter set for the generic model by Umlauf and Burchard (2003), where κ = 0.4,
− is based on the Algebraic Reynolds-Stress(cµ0 ) 2 = 0.3, and cψ1 = m are assumed. The value of the parameter cψ3
Model (ASM) of Canuto et al. (2001). In unstable situations (Pb > 0), a diﬀerent value of the parameter cψ3 needs to+ = 1.0 is adopted.
be used. In MRI.COM, cψ3

α
−2.0

L
0.20

m
1.00

ψ

σk
0.80

n
−0.67

σψ
1.07

cψ1
1.00

cψ2
1.22

−
cψ3
0.05

15.4.4 Stability functions
The version A of the Canuto et al. (2001) stability functions is expressed as follows:
0.1070 + 0.01741α N − 0.00012α M,
1 + 0.2555α N + 0.02872α M + 0.008677α 2N + 0.005222α N α M − 0.0000337α 2M
0.1120 + 0.004519α N + 0.00088α M
ĉµ′ = cµ′ (cµ0 ) 3 =,1 + 0.2555α N + 0.02872α M + 0.008677α 2N + 0.005222α N α M − 0.0000337α 2Mĉµ = cµ (cµ0 ) 3 =

where
αM =

E2 2
M ,
ϵ2

E2 2
N
ϵ2

αN =

(15.99)(15.100)(15.101)

are nondimensional parameters referred to as the shear number and the buoyancy number. Using (15.79), they are related
to the nondimensional parameter used by the Mellor and Yamada model (15.40) and (15.41) asgM =(cµ0 ) 6
2

α M,

gH =

(cµ0 ) 62

αN .

(15.102)

These are presented in Burchard and Bolding (2001)‡ . MRI.COM adopts these as the stability functions for the Umlauf
and Burchard (2003) model.

15.4.5 Implementation
Equations (15.83) and (15.84) are solved for the prognostic variable E and ψ in the subroutine gls_main in gls.F90 as
the last procedure of each time step. The forward finite diﬀerence is used in the time evolution. The implicit method isused for the vertical diﬀusion of the turbulent kinetic energy or generic variable and energy dissipation term, since theseterms could become significantly large (see Section 19.5). The new E and ψ are used to determine the coeﬃcients ofviscosity and diﬀusivity for the next time step.
The turbulent kinetic energy and the master length scale are defined at the center of the tracer cell (i, j, k − 21 ). The
specific expression for the discretized form of the turbulent kinetic energy (E) equation is as follows:
E n+11 − E n 1
k− 2

k− 2

∆t

1 [
=
∆zk− 1
2

K E k−1 (E n+13 − E n+11 )
k− 2

k− 2

∆zk−1

−

K E k (E n+11 − E n+11 ) ]
k− 2

1

n+1
n
n
+ Ps + Pb − (cµ0 ) 3 Ek−) 2 /l k−
1 (E
1.
k− 1
2

2

k+ 2

∆zk
(15.103)2′′′µ and cµ in Burchard and Bolding (2001) is diﬀerent than that of this chapter. ĉµ and ĉµ here correspond to cµ and cµof Burchard and Bolding (2001).

‡ Note that definition of c

– 186 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 15

Turbulence Closure Models

The specific expression for the discretized form of the generic variable (ψ) equation is as follows:ψ n+11 − ψ n 1
k− 2

k− 2

∆t

1 [
=
∆zk− 1Kψ k−1 (ψ n+13 − ψ n+11 )
k− 2

k− 2

−

∆zk−1

2

Kψ k (ψ n+11 − ψ n+11 ) ]
k− 2

k+ 2

∆zk

ψn 1
k−
1
n+1
n
n
2
+ n 2 (cψ1 Ps + cψ3 Pb ) − cψ2 (cµ0 ) 3 ψk−.
1 (E1 ) /l
k−
k− 21
E 1
2
2
k− 2

(15.104)

The shear and buoyancy production terms Ps and Pb are discretized as follows:1
Ps = K M k−1
2(un+13 − un+11 )(un+13 − un 1 )
k− 2

k− 2

k− 2

k− 2

∆zk−1 ∆zk− 1
2

1+ K M k−1
2
1
Pb = − K H k−1
2

1+ KM k2(v n+13 − v n+11 )(v n+13 − v n 1 )
k− 2

k− 2

k− 2

k− 2

∆zk−1 ∆zk− 1
2

B n+13 − B n+11
k− 2

k− 2

∆zk−1

1
− KH k2(un+11 − un+11 )(un 1 − un+11 )
k− 2

k+ 2

k− 2

∆zk

k− 2

k+ 2

∆zk ∆zk− 1
2

1+ KM k2

B n+11 − B n+11

k+ 2

(v n+11 − v n+11 )(v n 1 − v n+11 )
k− 2

k+ 2

k− 2

k+ 2

∆zk ∆zk− 1,(15.105)2.

(15.106)

The discrete expression for the shear production and the buoyancy production follows Burchard (2002), which is consistent
with the conservation law of the sum of mean and turbulent energy.
To summarize, the numerical operations proceed in the following order:
1. Update the turbulent kinetic energy and generic variable using (15.103), (15.104), (15.105) and (15.106) under theboundary conditions, (15.93), (15.94), (15.97) and (15.98).
2. Diagnose the turbulent length scale using (15.78),
3. Calculate the vertical viscosity and diﬀusivity using (15.88) and (15.88) for the next time step.

15.4.6 Limiters
Some kind of limiters for the turbulent variables may be needed to ensure numerical stability. MRI.COM supposes thefollowing conditions.
• Lower limiters of the turbulent kinetic energy and generic variable, Emin, ψmin , are introduced as basic limiters.
• If one of E and ψ becomes less than the limiter, turbulence is considered as stagnant. That is,
E = Emin,ψ = ψmin,l = l stag,

if either E ≤ Emin or ψ ≤ ψmin(15.107)• When turbulence is stagnant, the diagnosed viscosity and diﬀusivity (K M , K H ) should not exceed the backgroundvalues (K M bg, K H bg = 1.0 × 10−6 m2 s−1 ):K H = cµ′ (2Emin ) 1/2 l stag ≤ K H bg

K M = cµ (2Emin ) 1/2 l stag ≤ K M bg,(15.108)• Stagnant value of the turbulent length-scale (l stag ) is supposed to be the local vertical grid size (∆zlocal ) for numericalstability.
Table15.6

Default limiters for GLS model in MRI.COM.

1
2

(

Emin
10−6
0.1lstag)2

ψminl stagmin(cµ0 ) 2 E2/3

∆zlocallstag

15.4.7 Usage
Umlauf and Burchard (2003) model is invoked by GLS option on compilation.
Initial state for vertical diﬀusion of turbulence kinetic energy (K E ) and generic variable KΨ , turbulence kinetic energy(E), generic variable (Ψ), and length scale (l) is either read from restart file or set as the stagnant state defined in the– 187 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201715.4 Generic length scale model by Umlauf and Burchard (2003)
previous subsection. In the stagnant state, both of K E and KΨ are set as 1 × 10−6 m2 s. The initial state is specified bynamelist nml_gls_run. Parameters are listed on Table 15.7.
Tuning parameters may be specified at run time by namelist nml_gls. Parameters are listed on Table 15.8.
Table15.7 namelist nml_gls_run

variable namel_rst_gls_inunits
logicaldescription.true. : Read restart files specified by
nmlrs_gls_(avk, avq, eke, psi, alo) for the initial condition.
.false.: Start condition is that of no turbulence. See text for
details.

usage
default = l_rst_in of
nml_run_ini_state

Table15.8 namelist nml_gls

variable namemfnk
l_wave_damp_seaiceunits
1logicaldescription
η of Eq. (15.93)
Eq. (15.93) is used instead of (15.95)– 188 –

usage
default = 100.0
default = .false.

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 16

Bottom Boundary Layer (BBL)
MRI.COM has an option of adopting a simple bottom boundary layer model (BBL option). This chapter describes theformulation and usage of this model.

16.1 General description
In general, neither dense overflows from the Nordic-Scotland-Greenland ridges to the Atlantic Ocean nor the dense
downslope flows from the continental slope around Antarctica can be properly reproduced by z-coordinate models. The
former becomes the core water of the North Atlantic Deep Water, and the latter becomes that of the Antarctic Bottom
Water. As a result, the abyssal waters in the modeled ocean tend to exhibit a warming bias. To mitigate these deficiencies,
MRI.COM incorporates a simple bottom boundary layer (BBL) model as introduced by Nakano and Suginohara (2002).
In this model, BBL is attached immediately below the bottom of the interior oceanic domain and conceptually treatedas a plane along the sea floor at the model’s lowest vertical level (Figure 16.1). The following processes relevant to theoverflow/downslope flows are (partly) incorporated into this BBL model:• The advection along the bottom topography.
• The pressure gradient terms when the dense water lies on slope.
• The eddy activity to create the flow crossing f /h contours.

16.2 Grid arrangement and exchanges
In z-coordinate models, the flow along the bottom topography is expressed as a sequence of horizontal and verticalmovements along a staircase-like topography. When the number of vertical grid points representing the bottom topography
is roughly the same as that of the horizontal grid points, the downslope flow can be reasonably represented even in zcoordinate models (Winton et al., 1998). In general, however, the number of the vertical grids used in ocean generalcirculation models is not large enough. Even if the number is large enough, the concentrations of tracers might besignificantly diﬀused during their movement owing to the horizontal mixing with the surrounding waters. The incorporationof BBL remedies these problems.
The BBL grid cells in MRI.COM are arranged as in Figure 16.1. General interior grid cells in MRI.COM are composedof U-cells, but it is intuitively easier to understand the grid arrangement of BBL as if they are located at T-cells (represented
by "•" points of Figure 16.1a and b). The T-cells in BBL are located just below the bottom of interior T-cells even if aninterior T-cell consists of only one quarter of a U-cell (Figure 16.1b). The U-cells in BBL are regarded as slanting cells
connecting neighboring T-points in BBL (Figure 16.1d).
For the purpose of calculating horizontal exchanges between BBL cells due to advection and subgrid-scale mixing, theprimary BBL is attached below the normal ocean grid (Figures 16.1a and c). Assuming that the number of the verticalgrids, km, is 50 before BBL option is applied, the number of vertical grids, km, becomes 51 (50 + 1) when BBL option isapplied. In addition, to easily express the vertical exchange between the BBL cells and the inner ocean cells, we place adummy BBL cell just below the bottom grid of each inner column and copy the temperature, salinity, and velocities fromthe primary BBL cell (Figures 16.1b and d).
Regarding vertical advection, primary vertical mass transport at the upper surface of BBL is defined at T-cell as in theinterior (there is not vertical transport at the bottom of BBL). Even when the horizontal area of the bottom of the interior
T-cell is smaller than that of BBL, all vertical mass transport at the top of BBL flows into or out of the bottom of theinterior cell (Figure 16.1b and d). Tracers immediately above and below are carried by this vertical transport. In the samemanner as explained in Chapter 5, the continuity equation for the U-cell is defined as an average of those of surrounding– 189 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201716.3 Pressure gradient terms
T-cells (5.21). At each of the T-points connecting the bottom of the interior and the BBL, velocity (momentum) at all the
U-cells in BBL that surround the T-point participates in the vertical momentum advection.
Regarding vertical mixing, tracer in BBL is exchanged with the one at the bottom of the interior as in vertical advection(the blue arrow in Figure 16.1b). On the other hand, velocity in BBL is exchanged with the one at the bottom of the watercolumn of the same U-point (the blue arrow in Figure 16.1d), whose quarter cells may not necessarily constitute bottom
T-cells that communicate with T-cells of the BBL.

(a)

(b)

(w*,kv)(u*,v*)(w*,kv)(u*,v*)
k=km

k=km

(c)

(d)

(w*,kv)(w*,kv)

k=km

k=km

Figure16.1 Grid arrangement and schematic of horizontal and vertical exchange of tracers. (a, b) Arrangement ofT-cell represented by "•" symbols. (c, d) Arrangement of U-cell represented by "×" symbols. The primary BBL cells
lies at the bottom of the main (interior) oceanic domain (a and c). Horizontal advection and subgrid-scale mixing iscomputed there (solid black arrows). In relation to the interior, the T-cells in BBL are located just below the bottom ofinterior T-cells even when an interior T-cell consists of only one quarter of a U-cell (Figure 16.1b). The U-cells in BBL
are regarded as slanting cells connecting neighboring T-points in BBL (dashed blue line in d). Vertical transport intoand out of BBL occurs according to the primary vertical velocity defined at T-points (solid red arrows). Diagonally
upward / downward advection of momentum may occur along the dashed red arrows in (d). Vertical mixing of tracers
and momentum occurs between the cells indicated blue arrow in (b) and (d), respectively.

16.3 Pressure gradient terms
a. Formulation

BBL of MRI.COM is attached at the bottom of interior T-cells. U-cells of BBL are regarded as slanting and connecting
the surrounding T-points. (Figure 16.1d). In other words, sea floor is considered to have gradient along U-cells.
The pressure perturbation term in eq. (2.56) is1
1
1∇s p ′ =∇p′ s=−H (x,y) −∂s p′∇H (x, y)
ρ0
ρ0
ρ0
1
g ρ′
∇p′ s=−H (x,y) +z s ∇H (x, y),=ρ0ρ0
and the geopotential term is

g ρ′
g ρ′
g ρ′
∇s z =∇z s=−H (x,y) −z s ∇H (x, y).
ρ0ρ0ρ0– 190 –

(16.1)(16.2)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 16

Bottom Boundary Layer (BBL)

Then, the pressure gradient for a U-cell in BBL is expressed as
1
g ρ′
1
g ρ′
∇s p ′ +∇s z =∇p′ s=−HBBL (x,y) +∇z s=−HBBL (x,y) .
ρ0
ρ0
ρ0
ρ0

(16.3)

b. Discretization

Horizontal gradient in (16.3) is evaluated along the bottom slope, not along the constant s-surface.
The depth of U-point in BBL (z = zBBL ) that appears in the geopotential term is separated into time-independent(−HBBL ) and time-dependent (z ′) part, so that
g ρ′
g ρ′
g ρ′ ′∇z s=−HBBL (x,y) = −∇HBBL +
∇z s=−HBBL (x,y) .
ρ0
ρ0
ρ0

(16.4)The pressure perturbation at T-point in BBL, whose vertical grid index is k = ktbtm, is computed by adding an increment
due the downward vertical integration along the T-point from the corner of the upper face of the U-cell at the bottom ofinterior, whose vertical grid index is k = kbtm − 1, to the T-cell in BBL (k = ktbtm). The finite diﬀerence form is given as
follows:
ktbtm−1∑(dzt)i, j,ktbtm+ 1
2∆pi,′ j,ktbtm =( ρ′ )i, j,l+ 1 (dzt)i, j,l+ 1 + g( ρ′ )i, j,ktbtm+ 1.
(16.5)
2222
l=kbtm−1
The gradient of pressure perturbation is computed as′′(∇pBBL )i+ 1 , j+ 1 =(∇r p )i+ 1 , j+ 1 ,kbtm−1 + x̂
2

2

2

′
′
∆pi+1,∆p′+∆p′
j+1,ktbtm +∆pi+1, j,ktbtm− i, j+1,ktbtm2 i, j,ktbtm
2∆x i+ 1 , j+ 1

2

2

+ ŷ

2

′
′
∆p′
∆pi+1,+∆p′
j+1,ktbtm +∆pi, j+1,ktbtm− i+1, j,ktbtm2 i, j,ktbtm
2∆yi+ 1 , j+ 1
2.

(16.6)

2

The geopotential term is computed as[ g ρ′
ρ0

∇zBBL

=−

g
ρ0

]

ρ′i+1, j+1 +ρ′i, j+1 HBBL i+1, j+1 −HBBL i, j+1
ρ′i+1, j +ρ′i, j HBBL i+1, j −HBBL i, j+2∆xi+ 1 , j
2∆xi+ 1 , j+1
{[
2

g
ρ0

2

2
+

+

i+ 21 , j+ 12x̂ρ′i, j+1 +ρ′i, j HBBL i, j+1 −HBBL i, j
ρ′i+1, j+1 +ρ′i+1, j HBBL i+1, j+1 −HBBL i+1, j+2∆y2∆yi+1, j+ 1[
i, j+ 1
2

2

2

′
′
′
′
ρ′i+1, j +ρ′i, j zi+1,
ρ′i+1, j+1 +ρ′i, j+1 zi+1,
j −z i, j
j+1 −zi, j+1+2∆xi+ 1 , j
2∆xi+ 1 , j+1
{[
2

2

2

]
x̂ +

]

ŷ

]}

′
′
′
′
ρ′i, j+1 +ρ′i, j zi,
ρ′i+1, j+1 +ρ′i+1, j zi+1,
j+1 −z i, j
j+1 −z i+1, j+2∆yi, j+ 1
2∆yi+1, j+ 1
[
2

2

2

]}
ŷ .
(16.7)In the source code, the depth of T-point hot is used for HBBL . To be exact, (hot − dzt_cnst/2) should be used.
However, because the initial thickness of the BBL (dzt_cnst) is spatially uniform in MRI.COM, this omission would notaﬀect the result.

16.4 Eddy eﬀects
In the rotational frame, the dense water along the continental slope flows along f /h without the eddy eﬀects. In this case,introducing the BBL model does not lead to a better representation of the overflow/downslope-flow in the Nordic Seas oron the continental shelf around Antarctica.
In the real world, eddy eﬀects create the cross f /h flow, resulting in the overflow/downslope-flows. Jiang and Garwood(1996) investigated the three-dimensional features of downslope flows using an eddy-resolving model and demonstrated
that the dense water descends roughly 45° left to the geostrophic contour (in the northern hemisphere) with vigorous eddy
activity. The observation of significant eddy activity south of the Denmark Strait is consistent with this result. In order to– 191 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201716.5 Usage
incorporate this eﬀect into the non-eddy resolving models, we apply Rayleigh drag whose coeﬃcient is nearly equal to the
Coriolis parameter, α ≃ f . In this case, the geostrophic balance is modified and written as follows:
1 ∂p
− αu,
ρ0 ∂ x1 ∂p
fu = −− αv,
ρ0 ∂ y−fv = −

(16.8)(16.9)

After some algebra we obtain
(
)
α
1 ∂p
− f v − u = − f v′ = −,fρ0 ∂ x)(1 ∂pα,
f u + v = f u′ = −fρ0 ∂ y(16.10)(16.11)

where v ′ = v − (α/ f )u and u ′ = u + (α/ f )v correspond to the geostrophic velocity for the pressure gradient. If we put
α = f , the direction of the flow is 45° to the right of the pressure gradient. In general, the horizontal pressure field isparallel to the topographic contour near the deep water formation area. Thus, incorporating the Rayleigh drag coeﬃcient,
α = f , leads to the dense water descending 45° to the left of the geostrophic contour.
Because this Rayleigh drag is thought to be caused by the eddy activity and is observed where the dense water descends
to the deep layer, the coeﬃcient of the Rayleigh drag should be parameterized as a function of the local velocity andtopography. However it is very diﬃcult to appropriately determine this function in coarse-resolution models. Accordingly,
in the default setting of BBL option in MRI.COM, the depth range where α = f is arbitrarily set above 2500 m in thenorthern Atlantic and above 4000 m around Antarctica to represent the observed dense overflow/downslope-flow. Below
those depths, α is set to zero. This setting may be explicitly specified by using namelist nml_bblrayfrc (Table 16.1).

16.5 Usage
In this section, we show how to use the BBL model in the MRI.COM.
• Add BBL option in the file configure.in.
• Set kbbl = 1 in the file configure.in. (kbbl is the number of BBL layers. Presently, only one layer (kbbl = 1)
is supported).
• Set the lowest layer of the array dz(km) as the thickness of BBL. Fifty to one hundred meters is recommended forthe thickness of BBL, which roughly corresponds to the observed thickness of the BBL layer near the Denmark
Strait.
• Append ho4bbl, exnnbbl to file_topo after ho4, exnn (Section 21.2.2). These two variables have the sameformat as ho4, exnn. Variable ho4bbl corresponds to the thickness of the BBL, and variable exnnbbl corresponds
to the number of BBL layers; only one layer is presently supported. Thus exnnbbl should be equal to or less thanone.
• Set the region where non-zero Rayleigh drag coeﬃcient is applied (nml_bblrayfrc; Table 16.1). The default
setting is above 2500 m in the North Atlantic and above 4000 m around the Southern Ocean.
If model integration starts from the 3-D distribution of temperature and salinity, the temperature and salinity in the BBL
cells are set to those in the lowermost inner cells.

Example of configure.in with BBL option
OPTIONS="BBL"
IMUT=184
JMUT=171
KM=48
KSGM=48
KBBL=1



The file that contains topographic information is read by the model at run time as follows (see also Section 21.2.2).

– 192 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 16

Bottom Boundary Layer (BBL)

How topography data (file_topo) with BBL is read by model
integer(4) :: ho4(imut,jmut),exnn(imut,jmut)
integer(4) :: ho4bbl(imut,jmut),exnnbbl(imut,jmut)
integer(4),parameter :: lun = 10
open(lun, file=file_topo, form=unformatted)
read(lun) ho4, exnnread(lun) ho4bbl, exnnbbl ! This is for BBL.



Table16.1 Namelist nml_bblrayfrc for setting the lower limit above which Rayleigh damping is applied in BBL

variable
deprfcnunits
cmdeprfcscm

description
subjective lower limit of depth (for northern hemisphere) above which Rayleigh damping is applied
subjective lower limit of depth (for southern hemisphere) above which Rayleigh damping is applied

usage
default = 2500 m= 2500 × 102 cm
default = 4000 m= 4000 × 102 cm

16.6 Usage notes
16.6.1 Limit of the region where BBL model should be applied
Ideally, the BBL model should be applied universally and its behavior should change according to the local oceanic state.
In the real configuration, the BBL model is only eﬀective near the abyssal water formation areas. When used in otherareas, such as the near the equator, the BBL model does not improve the tracer and velocity fields. Furthermore, because
the BBL model connects the model cells along the topography, it inevitably induces a water mass exchange between thecold abyssal water and the warm shallow water. Thus, unphysical diapycnal diﬀusion could occur with the BBL model.
This eﬀect is not severe in high latitudes where the diﬀerence in temperature between the shallow continental shelves andthe deep layer is expected to be small, but it is extremely problematic for the cells in low latitudes. This problem is similar
to the problem for typical σ-layer models. To prevent this, it is recommended that the BBL model should be applied only
to the regions where the processes reproduced by the BBL model are thought to be important.

16.6.2 Limits of the BBL
Linear interpolation of the temperature and salinity along an extremely steep slope may cause problems. For the default
setting of MRI.COM, the BBL model is not applied in such places and isolated grids. This setting may be modified byusing namelist nml_bblgeo_limit (Table 16.2).
Table16.2 Namelist nml_bblgeo_limit for limiting the range where BBL scheme is applied

variable
hbblminunits
cm

bbldmaxcm

description
for sea floor depth shallower than hbbldmin, BBL
is not applied
if the sea floor is deeper than that of the neighboring point by bbldmax, the BBL is not applied

usage
default = 0 m
default = 10 000 m
10 000 × 102 cm

16.6.3 Notes for program coding and analyses
In each vertical column, the BBL exists both at the bottom cell (k = km) and at the lowest ocean cell (k = exn(i, j)). In
general, the lowest ocean cell of the T-point and U-point for the same horizontal indices might diﬀer. Thus, the treatmentof topographic mask array atexl and aexl is always very confusing. At the bottom cells (k = km) that corresponds theprimary BBL cells, we set atexl = 1 and aexl = 1 while we set atexl = 1 and aexl = 0 at the lowest ocean cells that
corresponds the dummy BBL cells. In taking global average, the treatment of atexl = 1 in the BBL model needs special
care to avoid double counting.

– 193 –

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 17

Sea iceThis chapter describes the sea ice part. The sea ice part of MRI.COM treats formation, accretion, melting, and transfer
of sea ice and snow. Heat, water, salt, and momentum fluxes are exchanged with the ocean. Sea ice is categorized by its
thickness, but it has a single layer. Snow does not have heat capacity (so-called zero-layer). Thus, it might be regarded as
an intermediate complexity ice model.
This chapter is organized as follows. Section 17.1 outlines the model. The following sections describe details of
the solution procedure. According to the order of solving the equations, we deal with thermodynamics in Section 17.2,
remapping among thickness categories in Section 17.3, dynamics in Section 17.4, advection in Section 17.5, and ridging
in Section 17.6. Adjustment on the sea ice distribution is explained in Section 17.7. Discretization issues are described in
Section 17.8 and some technical issues are presented in Section 17.9. Finally, usage notes are presented in Section 17.10.

17.1 Outline
The sea ice part of an ocean-sea ice model gives sea surface boundary conditions to the oceanic part. Heat, fresh water,
salt, and momentum are exchanged at their interfaces. The sea ice part solves fractional area, heat content, thickness, andtheir transport of ice categorized by its thickness and dynamics of the grid-cell averaged ice pack.
The ice model of MRI.COM is based on the ice-ocean coupled model of Mellor and Kantha (1989). For processes thatare not explicitly discussed or included there, such as categorizing by thickness, ridging, and rheology, we adopt those of
the Los Alamos sea ice model (CICE) version 3.14 (Hunke and Lipscomb, 2006).
The fundamental property that defines the state of sea ice is the fractional area as a function of location (µ, ψ) andthickness (h I ). The equation for this distribution function (g(µ, ψ, h I )) is∂g∂
1 ( ∂(ghψ u I ) ∂(ghµ vI ) )=−( f g) −
+
+ χ,
∂t
∂h I
h µ hψ
∂µ∂ψ(17.1)

where f is the thermodynamic growth rate of thickness, (u I , vI ) is the velocity vector of ice pack, and χ is the rate ofchange of distribution function caused by ridging. The growth rate of ice thickness is computed by solving thermodynamic
processes. Using this rate ( f ), thickness categories are remapped according to the first term on the r.h.s. of (17.1). To
compute the velocity of the ice pack (u I , vI ), we have to solve the momentum equation. On transporting the ice distribution(second and third terms on the r.h.s.), other conservative properties such as volume and energy are also transported. Using
the transported ice distribution function, the ridging process ( χ; fourth term on the r.h.s.) is solved. The formulation andsolving procedure of each process are presented in later sections.
We discretize the thickness in several categories. If an ice pack is divided into n categories separated at Hn with
H0 = 0 m, the fractional area of each category (an ) is defined as follows:
an =

∫ Hn
gdh.

(17.2)

Hn−1

Other major variables, ice and snow thickness, surface temperature, bottom temperature and salinity, and internal energyof ice, are defined for each category. Velocity is defined for an ice pack, the total ice mass in a grid cell. In the verticaldirection, both sea ice and snow have one layer. Sea ice has heat capacity, but snow does not. The heat capacity for sea ice
is due to brine and is represented by the temperature at the center of the sea ice. Figure 17.1 and Table 17.1 summarize
symbols used in this chapter and their variable names in the source code of MRI.COM. Note that the symbol T is used torepresent the temperature of any substance treated in this chapter. The unit is ◦C, though the use of an upper case letter(T) may recall absolute temperature.

– 195 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.1 Outline
Table17.1 Physical quantities used in the sea ice part (cf. Figure 17.1) and their variable names in the source code.

hI
hs
A
Ah I
Ahs
E1
Ah I E1
T3
T2
T1
T0 I (= T0 )T0 L
S0 I (= S0 )
S0 L
QI2
QIO
QAO
FTI
FTL
FS I
FS L
W
WAI
WAI
WIO
WAO
WROiceMeaning

Units

Array name foreach thickness
category

Array name foraverage or sum
over thickness
categories

ice thickness
snow thickness
area fraction (compactness)
grid cell average of ice thickness (volume per unit area)
grid cell average of snow thickness (volume per unit
area)
enthalpy of ice
grid cell average of enthalpy
skin temperature of upper surface
temperature at snow-ice interface
temperature of ice
skin temperature of lower surface
skin temperature of sea surface at open leadsskin salinity of lower surface
skin salinity of sea surface at open leads
heat flux in the upper half of ice (= Q S )heat flux on the ice side of the ice bottom
heat flux on the air side at open leads
heat flux on the ocean side of the ice bottom
heat flux on the ocean side at open leads
virtual salinity flux below sea ice driving the top layerof the ocean modelvirtual salinity flux in open water driving the top layerof the ocean model
fresh water flux driving the top layer of the ocean model(= AWIO + (1 − A)WAO + WRO + WFR )fresh water flux due to snow fall at the upper surface ofice
fresh water flux due to sublimation at the upper surfaceof icefresh water flux due to freezing and melting at the bottomof icefresh water flux due to freezing at open leads
fresh water flux due to melting of sea ice at the uppersurface of ice

m
m
1
m
m

hicen
hsnwn
aicen
hin
hsn

hiceo
hsnwo
a0iceohihsnw

J kg−1J kg−1 m
◦C
◦C
◦C
◦C
◦C
pss
pss
W m−2
W m−2
W m−2
W m−2
W m−2
pss m s−1

qicen
eicen
tsfcint2icent1icent0icen(1 :)t0icen(0)s0n(1 :)s0n(0)
fheatufheatn
fheat
ftio
ftao
-

tsfci
t0iceot0icel-pss m s−1

-

-

m s−1

-

wfluxim s−1-

snowfallm s−1

sublim

-

m s−1

wio

-

m s−1
m s−1

wao
-

wrsi

Continued on next page

– 196 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea ice
Table 17.1 – continued from previous page
Meaning
Units

WROsnow

fresh water flux due to melting of snow at the uppersurface of ice
fresh water flux due to formation of frazil ice
zonal component of ice pack velocity
meridional component of ice pack velocity

WFR
uI
vI

Array name foreach thickness
category

Array name foraverage or sum
over thickness
categories

m s−1

-

wrss

m s−1
m s−1
m s−1

-

wrso
uice
vice

shortwave

QAI

WAI
T3

snow

hS

QS
T2

Ss = 0.0

S0L
FSL

FTL

QI2

ice

WAO

QAO
T0L

E1

SI = 4.0

T1

hI
QIO

S0(=S0I)

T0 (= T0I)
FTI

FSI

T(1), S(1)

WFR

WRO
WIO

Figure17.1 Meaning of symbols and their locations. The left side is related to heat flux, and the right side is related
to fresh water flux. A column of ice is separated into the part that originated from sea water (thickness: h I and salinity:
SI = 4.0 pss) and the part that originated from snow (thickness: hs ). The former is further divided into upper andlower halves. Thus a column of ice has three vertical layers. The temperatures at the layer boundaries are T0 , T1 , T2 ,and T3 from the bottom. Heat fluxes within each layer are QIO , QI2 , and Q S from the bottom. The heat flux at the
air-ice interface is QAI and that at the ice-ocean interface is FTI . Sea ice is in fact categorized by thickness, and eachsymbol should have a suﬃx (n) of the category number. At open leads or open water, symbols have the suﬃx L. For
the definitions of fresh water fluxes, see Table 17.1.

17.2 Thermodynamic processes
Although the surface fluxes are positive downward (positive toward the ocean) in the ocean model, the sea ice part is codedsuch that fluxes are positive upward. In this section, we assume that the fluxes are positive upward.

17.2.1 Thermal energy of sea ice
In considering thermodynamics, thermal energy of sea ice should be defined. The basis of energy (i.e., zero energy) isdefined here as that of sea water at 0 ◦C. The thermal energy (enthalpy; E(T, r)) of sea ice that has temperature T (< 0)[◦C]and brine (salt water) fraction r is the negative of the energy needed to raise the temperature to 0 ◦C and melt all of it:
E(T, r) = r (CpoT ) + (1 − r)(−L F + CpiT ),(17.3)

where Cpo and Cpi are the specific heats of sea water and sea ice, and L F is the latent heat of melting/freezing. Thus
defined, the energy of sea ice is negative definite. The positive brine fraction represents the fact that the salt with salinity(SI ) in sea ice exists in a liquid state. If the ice temperature is T1 , the brine is assumed to have the same temperature, andits salinity is S = T1 /m, where m defines the freezing temperature as a function of salinity. Hence, the brine fraction ofsea ice is r = SI /S = mSI /T1 . As in Mellor and Kantha (1989), the specific heat of snow is not considered.
– 197 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.2 Thermodynamic processes
Once sea ice enthalpy E is known, sea ice temperature T1 is determined by solving the quadratic equation,
E=

(
mSI
mSI )
(CpoT1 ) + 1 −(−L F + CpiT1 ),T1
T1

which gives√−b − b2 − 4ac
,
T1 =
2a

(17.4)(17.5)

where
a = Cpi, b = −L F − E + mSI (Cpo − Cpi ), c = mSI L F .

(17.6)

Figure 17.2 shows that the other root is unphysical.

Figure17.2 Sea ice enthalpy E (17.4) as a function of sea ice temperature (T1 ). Units in J kg−1 .

17.2.2 Formation of new sea ice
a. From sea water with temperature below the freezing point

Sea ice is formed when sea surface temperature is below the freezing point. If the temperature of the first layer of theocean model is below the freezing point as a function of salinity, the temperature is set to the freezing point (Tfreeze ), andthe heat needed to raise the temperature is regarded as the release of latent heat and is used to form new ice.
The thickness of the new ice (h I ) is computed by assuming that the total thermal energy of the first layer of the ocean
model (whose layer thickness is ∆z 1 ) is conserved before and after the sea ice formation:2

ρo ∆z 1 [CpoT] = ρ I h I {r (CpoTfreeze ) + (1 − r)(−L F + CpiTfreeze )} + ρo (∆z 1 − ρ I h I /ρo )(CpoTfreeze ),22(17.7)

where r = mSI /Tfreeze . Using the above equation, we compute the thickness of the new ice:hI =ρo ∆z 1

Cpo (Tfreeze − T ).
ρ I (1 − r) L F + Tfreeze (Cpo − Cpi )2(17.8)

With the ice thickness known, the fractional area is determined by the following procedure:• For a grid cell without sea ice, the first thickness is set to hnew
= 10 cm and the fractional area is computed asI
new = h and A = 1.
A = h I /hnew
.
If
A
>
1
(i.e.,
h
>
10
cm),
h
I
I
I
I
• For a grid cell where sea ice already exists, h I is added to each category and open water.
– 198 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea ice

Note that this operation practically eliminates super-cooling in the ocean interior. Hence, the formation of frazil ice is notconsidered in this model.
The medium that contacts the sea ice diﬀers at the upper (air) and lower (sea water) surfaces. Even for a grid cellwith sea ice, not all the surface is covered by sea ice, i.e., there may be open water. We treat processes at each interface
separately. Heat flux and thermal energy are computed by solving the balance equation at each interface and each layer inFigure 17.1.
b. Input of iceberg

Input of iceberg (Ficeberg ) is given as a water mass flux per unit area (kg m−2 s−1 ). The thickness of the new ice (h I ) is
given byh I = Ficeberg ∆t/ρ I ,(17.9)
where ∆t is the unit time step. With the ice thickness known, the fractional area is determined by the same procedure used
for ice created from supercooled water as explained above.

17.2.3 Air-ice interface
a. Heat flux at the upper surface of ice (Q AI )The surface heat flux (QAI ) is expressed as follows:
QAI = QSI + QLI − (1 − α I )SW − LW + ϵ I σ(T3 + 273.15) 4,(17.10)where QSI is sensible heat flux, QLI is latent heat flux, SW is downward shortwave radiation flux, LW is downward
longwave radiation flux, α I is albedo, ϵ I is emissivity, and σ is the Stefan-Boltzmann constant. We will now examine
each component.
The downward shortwave radiation is represented by SW. The albedo of sea ice is α I , which is 0.82 forcold (T3 < −1 ◦C) snow, 0.73 for melting snow, and 0.64 for bare ice (while melting).
One might use a more sophisticated albedo scheme included in the Los Alamos sea ice model (CICE; Hunke andLipscomb, 2006) by choosing CALALBSI option. We briefly describe how the incoming shortwave radiation is treated bythis scheme. The downward shortwave radiation is treated at each interface as follows:■ i. Shortwave:• Among the net absorbed shortwave flux (= (1 − α I )SW), some fraction (i 0 ) penetrates into ice and the rest isabsorbed at the surface and used to warm the upper interface. See Table 17.2 for the specific value of i 0 .
• The part penetrating into the ice (= (1 − α I )i 0 SW) is attenuated according to Beer’s Law with the bulk extinction
coeﬃcient κ i = 1.4 m−1 . The attenuated part is used to warm the ice interior.
• The rest enters the ocean.
The albedos and penetration coeﬃcients of CICE are listed on Table 17.2. The property f snow is the snow fraction of
the upper surface of the ice, which is expressed as follows:
f snow =hs,
hs + hsnowpatch(17.11)

where hs is the snow thickness and hsnowpatch = 0.02 m.
If the ice thickness (h I ) is less than href = 0.5 m, the albedo of thin ice is computed as
αthinice = α o + β(αcoldice − α o ),
where
β=

arctan(ar h I ),arctan(ar href )ar = 4.0,(17.12)(17.13)and α o is the albedo of the ocean.
If the surface temperature T3 becomes −1 < T3 < 0 ◦C, the albedo of melting ice and snow is computed as
αmeltice = αthinice − γi (T3 + 1.0),
αmeltsnow = αcoldsnow − γs (T3 + 1.0),– 199 –

(17.14)(17.15)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.2 Thermodynamic processes
where the condition αmeltice > α o is imposed. Using the snow fraction on the surface of the ice f snow , the total albedo iscomputed as
αi = αmeltice (1 − f snow ) + αmeltsnow f snow .
(17.16)
The albedos for visible and near infra-red wave lengths are computed separately. If the shortwave flux is given as thesum of all four components (direct and diﬀuse for visible and near infra-red wave lengths), a constant ratio (visible) : (nearinfra-red) = 0.575 : 0.425 is assumed, and the total albedo is computed as the weighted average.
Table17.2 Albedo and surface transparency of the albedo scheme of CICE.

albedo for cold snow αcoldsnow (T3 < −1 ◦C)
albedo for cold ice αcoldice (T3 < −1 ◦C, hi > 0.5 m)
reduction rate of albedo for melting ice γi (−1 ◦C < T3 < 0, hi < 0.5 m)
reduction rate of albedo for melting snow γs (−1 ◦C < T3 < 0, hi < 0.5 m)
fraction of transparent shortwave flux through the ice surface (i 0 )near infra-red(> 700 nm)
0.70
0.36
−0.075 /◦C−0.15 /◦C
0.0

visible(< 700 nm)
0.98
0.78
−0.075 /◦C−0.10 /◦C
0.7 × (1.0 − f snow )■ ii. Longwave: The downward longwave radiation from the atmosphere is represented by LW in (17.10). The black
body radiation from the ice surface is ϵ I σ(T3 + 273.15) 4 , where ϵ I is emissivity, and σ is the Stefan-Boltzmann constant.
Hereinafter, we use LW I = LW − ϵ I σ(T3 + 273.15) 4 as the net longwave radiation.
■ iii. Sensible heat flux:

The sensible heat flux (QSI in (17.10)) is computed using a bulk formula:
QSI = ρ a Cpa CHAI U10 (T3 − TA ),(17.17)

where ρ a is the density of air, Cpa is the specific heat of air, CHAI is the bulk transfer coeﬃcient for heat, U10 is the scalarwind speed at 10 m, TA is the surface air temperature, and T3 is the ice surface temperature (Figure 17.1).
■ iv. Latent heat flux:

The latent heat flux (QLI in (17.10)) is computed using a bulk formula:
QLI = ρ a L s CEAI U10 (qi − q A ),(17.18)where L s is the latent heat of sublimation, CEAI is the bulk transfer coeﬃcient for moisture, qi is the saturation humidity at
the ice surface temperature T3 , and q A is the specific humidity of air. Section 14.12.3 details a computing method for qi .
Fresh water loss due to sublimation is computed as
WAI = ρ a CEAI U10 (qi − q A )/ρo .

(17.19)■ v. Bulk coeﬃcient over sea ice:

Over sea ice, it is common to set constant neutral bulk coeﬃcients at 10 m, as describedin Large and Yeager (2004). Then, bulk coeﬃcients and wind speed at each height are calculated on the basis of stability.
Large and Yeager (2004) uses
CDAI10 = CHAI10 = CEAI10 = 1.63 × 10−3,(17.20)
but the default in MRI.COM uses, following Mellor and Kantha (1989),
CDAI10 = 3.0 × 10−3

(17.21)

CHAI10 = CEAI10 = 1.5 × 10−3 .

(17.22)

These coeﬃcients are given by namelist nml_air_ice at run time. See Table 17.9.
The wind stress over sea ice can be calculated with TAUBULK option.
b. Heat flux in the snow (Q S )If we neglect the heat capacity of snow, the heat flux is constant within the snow layer and is computed as
QS =ks(T2 − T3 ),hswhere hs is the thickness of the snow layer and k s is the thermal conductivity of snow.
– 200 –

(17.23)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea icec. Heat flux in the ice interior (Q I2, Q IO )In the upper half of the ice layer, the heat flux is computed as follows:
QI2 =

kI
(T1 − T2 ),h I /2(17.24)

where k I is the thermal conductivity of sea ice. If we neglect heat capacity of snow, Q S = QI2 . Using this relation, theinterface temperature is computed askskI
h s T3 + h I /2 T1
T2 = k
.
(17.25)kIsh s + h I /2
In the lower half of the ice layer, the heat flux is computed as follows:
QIO =

kI
(T0 − T1 ).
h I /2(17.26)d. Melting at the upper surface (WRO )The melting rate at the upper surface is computed as follows: First, the surface temperature (T3 ) is computed by equating
the fluxes at the upper surface (QAI = Q S ). If T3 is lower than the freezing temperature (0.0 ◦C for snow and mSI [◦C] forsea ice), melting does not occur. If T3 is higher than the freezing temperature, melting occurs. In this case, T3 is set to thefreezing temperature, and the interior heat flux just below the ice surface (Q S ) is recalculated. The imbalance between
QAI and Q S is used to melt snow or ice. The fresh water flux due to melting (by sign convention in this Chapter, this isalways negative) is computed as
WRO = (QAI − Q S )/( ρo L 3 )(17.27)
where L snow
= L F for snow melt and3L ice
3 ≡ [E(T3, 1) − E(T1, r 1 )](17.28)for ice melt. The brine fraction for sea ice is r 1 = mSI /T1 . The temperatures of the melted water are mSI [◦C] for sea ice
and 0 ◦C for snow.
If all pre-existing snow and ice is going to be melted, that is, (Q S − QAI )∆t > ρs hs L snow
+ ρ I h I L ice
3 , a new interior flux
3′
(Q S ) is computed as follows:′
Q S = QAI + ( ρs hs L snow
+ ρ I h I L ice(17.29)3 )/∆t.
3′

In this case, Q S is equated with QIO and treated as the sensible heat flux at the ice-ocean interface.
e. Notes• All precipitation on sea ice is assumed to be snow.
• Melted water is assumed to run oﬀ to the ocean immediately, that is, there is no melt pond.
f. Solution procedure

The solution procedures are basically the same with or without snow. To be exact, the interface fluxes should be computed
iteratively by adjusting surface temperature T3 until a balance is achieved. We adopt the semi-implicit method described
below. A situation without snow (hs = 0, T3 = T2 ) is considered. First, the surface temperature (T3 ) is computed byassuming that the fluxes on both sides are the same. Inserting (17.10) and (17.24) into QAI = QI2 with T3 → T3 + δT3 ,kI{T1 − (T3 + δT3 )} = QLI (T3 + δT3 ) + QSI (T3 + δT3 )h I /2− (1 − α I )(1 − i 0 )SW − LW + ϵ I σ{(T3 + δT3 ) + 273.15}4 .

(17.30)By expanding the specific heat, latent heat, and black body radiation in a Taylor series, we have,
∂QSI
kI
∂QLI
(T1 − (T3 + δT3 )) = QLI + QSI +δT3 +δT3 − (1 − α I )(1 − i 0 )SW
h I /2∂T3∂T3− LW I + 4ϵ I σ(T3 + 273.15) 3 δT3 .
– 201 –

(17.31)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.2 Thermodynamic processes
Using this, we compute δT3 and add it to T3 to obtain a new temperature:δT3 =−QSI − QLI + (1 − α I )(1 − i 0 )SW + LW I + hkI I/2 (T1 − T3 )
∂QSI
∂QLI
kI
3
∂T3 + ∂T3 + 4ϵ I σ(T3 + 273.15) + h I /2,

(17.32)andT3new = T3old + δT3,

(17.33)

where∂QSI= ρ a Cpa CHAI U10,
∂T3
∂qi
∂QLI
= ρ a L s CEAI U10
.
∂T3
∂T3

(17.34)(17.35)Note that the dependency of L s on temperature (Section 14.12.3) is not considered in the partial diﬀerentiation with respect
to temperature. The specific form for the partial derivative of specific humidity (∂qi /∂T3 ) is presented in Section 17.11.1.
If the new surface temperature (T3new ) is below the freezing point, melting does not occur. If not, T3new is set to thefreezing temperature (= mSI ), and the heat flux in the ice interior is re-evaluated. The amount of melting (WRO ) is obtainedusing the imbalance:
QAI = − QSI − QLI + (1 − α I )(1 − i 0 )SW + LW I ,kI(T1 − T3new ),QI2 =h I /2L 3 = [E(T3new, 1) − E(T1, r 1 )],(QAI − QI2 )∆t∆h I == WRO ∆t.
ρI L3

(17.36)

If there is snow, the above procedure is performed for the snow surface (T3new = 0 ◦C, L 3 = L F ). If all the snow melts
away, the residual heat (Eres ) is used to melt the ice:
Eres = (QI2 − QAI )∆t − hs ρs L F ,L 3 = [E(mSI , 1) − E(T1, r 1 )],∆h I = − Eres /ρ I L 3 = WRO ∆t.

(17.37)17.2.4 Heat balance in the ice interior
The thermal energy of the ice is aﬀected by vertical heat fluxes and horizontal heat transport due to advection. The
equation for the thermal energy (enthalpy) is written as follows:
ρI hI

[∂
∂t

E(T1, r 1 ) + u I i]∂
E(T1, r 1 ) = QIO − QI2 + [SW surface − SW bottom ],∂ xi(17.38)where
SW surface = (1 − α I )i 0 SW,
SW bottom = (1 − α I )i 0 SW × exp(−κ i h I ).

(17.39)(17.40)

The above equation can be solved explicitly without causing serious problems when the time step is not too long.
If the enthalpy exceeds the upper limit, E(mSI , 1), QIO is adjusted so that the enthalpy gets the upper limit value at theend of the time step (∆t),
′

QIO =

ρ I h I [E(mSI , 1) − E n ]− [SW surface − SW bottom ] + QI2,
∆t′

(17.41)

where E n is the enthalpy at the time level n. QIO is used instead of QIO as the sensible heat flux at the ice-ocean interface.

– 202 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea ice17.2.5 Ice-ocean interface
Melting and freezing at the ice-ocean interface is computed using heat fluxes at the interface as depicted in Figure 17.1.
The ice-covered area and the open water are treated separately. The solution method slightly diﬀers from that of Mellor
and Kantha (1989).
In the ice-covered area, the heat flux on the ice side of the ice-ocean interface (QIO ) is computed according to (17.26).
In open water, the heat flux on the air side of the air-ocean interface (QAO ) is computed as in Chapter 14,
4
Qall
AO = Q SO (T ) + Q LO (T ) − (1 − α o )SW − LW + ϵ o σ(T + 273.15) .

(17.42)The temperature at the first level of the ocean model (T) is used. All the heat and fresh water fluxes are evaluated using thetemperature and salinity at the first level of the ocean model. By doing so, the equation to compute melting and freezing
rates becomes linear.
Here, shortwave radiation is assumed to pass through the skin layer without absorption and is excluded from theevaluation of the freezing rate in open water:
QAO = QSO + QLO − LW + ϵ o σ(T + 273.15) 4 .

(17.43)This operation causes the shortwave radiation to be absorbed in the ocean interior. In reality, the heat stored in the skinlayer in open water is used to melt ice laterally (edge melting). To include this eﬀect, a certain fraction (Ψ) of the bottom
melting (WIO ) may be used for the edge melting. The details are described in the last part of this section.
Melting and freezing occur owing to the imbalance between fluxes above (QIO , QAO ) and below (FTI , FTL ) the interface:
FTI = QIO − WIO ρo L o,
FTL = QAO − WAO ρo L o,

(17.44)(17.45)L o ≡ [E(T, 1) − E(T1, r 1 )].

(17.46)

where
T is the temperature at the first level of the ocean model. The heat flux that drives the first level of the ocean model isgiven by
FT = ( AQIO + (1 − A)QAO ) − WO ρo L o,
(17.47)
where
WO ≡ AWIO + (1 − A)WAO .

(17.48)

Following the formulation adopted by Mellor and Kantha (1989), we introduce skin layer temperature and salinity (T0 I ,T0 L , S0 I , and S0 L ) to solve for WIO and WAO , and thus FTI and FTL . To incorporate skin layer salinity in the systemof equations, we consider the problem in terms of virtual salt flux, in which the eﬀect of freshwater flux on salinity is
considered in terms of salt flux by keeping the water volume. The flux balance for salt below the interface is written as
follows:
FS I = WIO (SI − S),
FS L = WAO (SI − S).

(17.49)(17.50)Here, unlike Mellor and Kantha (1989), the salinity at the first level of the ocean model (S) is used instead of the salinity
at the skin layer (S0 I , S0 L ). By doing so, the equations to solve for S0 I and S0 L become linear as shown below. It could
also be said that it is natural to use the first level salinity itself in evaluating the salt flux that drives the first level ofthe ocean model. Note that only fresh water fluxes that are relevant to freezing and melting at the ice-ocean interface
are included in the above equations. The restoration to climatological salinity and fresh water fluxes caused by surface
melting, precipitation, and evaporation are excluded in the above balance.
The virtual salinity flux caused by melting and freezing at the ocean surface that drives the first level of the ocean model
is given by
FS = ( AWIO + (1 − A)WAO )(SI − S).
(17.51)
Fluxes on the oceanic side of the interface (FTI , FTL , FS I , and FS L ) can also be obtained as the boundary conditions(z → 0) for the molecular boundary layer:
FTI /( ρo Cpo ) = − CTz (T0 I − T ),FTL /( ρo Cpo ) = − CTz (T0 L − T ),– 203 –

(17.52)(17.53)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.2 Thermodynamic processes
and
FS I = − CSz (S0 I − S),
FS L = − CSz (S0 L − S),
where

uτ
,
−1
(Prt k ln(−z/z0 ) + BT )CTz =(17.54)(17.55)(17.56)

2 + τ 2 ) 1/4 ρ−1/2 is the friction velocity, k = 0.4 is von Karman’s constant, z is the roughness parameter,uτ ≡ (τIO
0
o
IO y
x
(τIO x , τIOy ) is the stress vector at the ocean-ice interface, andBT = b( z0 uτ ) 1/2
ν

Pr 2/3,(17.57)

with Pr ≡ ν/α t = 12.9. The specific values for other parameters are given in Section 17.11.2.
Parameters related to salinity are given by
CSz =

uτ
,
(Prt k −1 ln(−z/z0 ) + BS )and
BS = b( z0 uτ ) 1/2
ν

Sc2/3,(17.58)(17.59)where Sc ≡ ν/α b = 2432.
Following Mellor and Kantha (1989), roughness parameter z0 is computed as follows:
lnz0 = Alnz0 I + (1 − A)lnz0 L ,where
z0 I = 0.05

hI,
h I lim(17.60)h Ilim = 3.0 m,(17.61)ρo uτ2.
ρa g(17.62)andz0 L = 0.016

However, the roughness parameter below ice (17.61) is also used for open water (z0 L ) in MRI.COM.
The above equations are solved simultaneously to obtain melting and freezing at the upper surface of the ocean under
the following constraints:
WO = 0, if A = 0,
(17.63)andT0 I = mS0 I

and T0 L = mS0 L , if A > 0,

(17.64)

where m defines the freezing line as a function of salinity.
Solution procedure

We first solve for S0 I and S0 L using (17.44), (17.49), (17.52), (17.54), (17.45), (17.50), (17.53), and (17.55):
CSz S + ( ρo Cpo CTz T − QIO )(SI − S)/ρo L o,CSz + ρ0 Cpo CTz m(SI − S)/ρo L o
CS S + ( ρo Cpo CTz T − QAO )(SI − S)/ρo L o
S0 L = z.
CSz + ρo Cpo CTz m(SI − S)/ρo L o
S0 I =(17.65)(17.66)Using S0 I and S0 L , T0 I and T0 L are computed from (17.64), FTI and FTL are computed from (17.52) and (17.53), FS I and
FS L are computed from (17.54) and (17.55), and finally WIO and WAO are computed from (17.44) and (17.45). Since WAO
should be positive by definition (only freezing is allowed in open water), WAO is set to zero when WAO < 0. In this case,
S0 L = S from (17.50) and (17.55). Furthermore, FTL = Q AO in (17.45) and T0 L = T assuming that there is no ice eﬀect inopen water if freezing does not occur. Finally, the surface boundary conditions for the ocean model, the surface heat andvirtual salt fluxes, are computed from (17.47) and (17.51). If we allow variation of volume of the first layer of the ocean– 204 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea icemodel, virtual salt flux is not necessary. The fresh water flux (17.48) is used as a boundary condition, assuming that WO
has salinity of SI and temperature of the first layer of the ocean model.
If −ρo WIO ∆t/ρ I > h I , all ice would melt away. In this case, the amount of heat needed to melt all the ice is consumed,
and the residual is returned to the ocean. That is, the heat flux necessary to melt all the ice is removed from the ocean.
′Specifically, WIO is recomputed byρo ′h I = − WIO ∆t(17.67)
ρI
and FTI is obtained as′

FTI = QIO − WIO ρo L o .

(17.68)

If snow remains above the completely melted sea ice, the snow is also melted using ocean’s heat, which is treated as thelatent heat loss from the ocean (QLO ),QLO = ρs hs L F .
(17.69)When −ρo WIO ∆t/ρ I < h I , reduction of the fractional area is allowed by edge melting. The procedure is as follows:
The fraction used for edge melting is defined as Ψ, and −(1 − Ψ)WIO ρo ∆t/ρ I is used to reduce the thickness holding the′fractional area. The new thickness is h I and Ψ ρo WIO ∆t/ρ I is used to reduce compactness holding the thickness:′

A′ = A(1 + Ψ ρo WIO ∆t/ρ I h I ).

(17.70)

There seems to be no widely accepted parameterization scheme for edge melting. According to Steele (1992), bottom
and top melting are dominant processes, and thus Ψ ∼ 0.0 − 0.1 is usually used in MRI.COM.

17.2.6 Archimedes’ Principle
From Archimedes’ principle, the part of snow that is below freeboard absorbs sea water to become sea ice. The following
equality will be achieved at equilibrium:ρ I h I + ρ s hs = ρo h I .
(17.71)At the end of the time step, the above equality is checked, and if
ρ I h I + ρs hs> hI ,
ρ0(17.72)hnew=
Iρ I h I + ρs hs.
ρo(17.73)δhs = −(hnew
− hI ) ρII,ρs(17.74)then the new ice thickness is set as:

Since the change in snow thickness isthe new snow thickness is obtained as follows:hnew= hs + δhs .
s(17.75)Since the salinity of sea ice is SI = 4.0 pss, salt is removed from the first layer of the ocean model:
S new = S old − SI

δhs ρs.
∆z1 ρo(17.76)17.3 Remapping in thickness space
After the thermodynamic processes are solved, the resultant ice thickness in some thickness categories might not be within
the specified bound. Following the method adopted by CICE, we assume that there is a thickness distribution function ineach category and use it to redistribute the new thickness distribution into original categories.
This procedure corresponds to the first term on the r.h.s. of (17.1):∂∂g( f g).
=−∂t∂h I– 205 –

(17.77)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.3 Remapping in thickness space
In practice, a thickness category is regarded as a Lagrange particle, and the category boundaries are displaced as aresult of thermodynamics. A linear thickness distribution function is assumed within each displaced category, and ice isremapped into the original categories using these functions.
First, boundaries of thickness categories are displaced. If the ice thickness in category n (hn ) changes from hnm to hnm+1(m is the time step index), the growth rate ( f n ) at thickness hn is represented as f n = (hnm+1 − hnm )/∆t. Using this, thegrowth rate (Fn ) at the upper category boundary Hn is obtained by linear interpolation:
Fn = f n +f n+1 − f n(Hn − hn ).
hn+1 − hn(17.78)

If the fractional area is zero in either category n or n + 1, Fn is set to the growth rate at the non-zero category. When thefractional area is zero on both categories, Fn = 0. The new category boundary after thermodynamics is obtained as
Hn∗ = Hn + Fn ∆t.

(17.79)∗ , H ∗ ] is determined. For simplicity, we
Next, the thickness distribution function g within the displaced category [Hn−1n∗∗write H L = Hn−1 and HR = Hn . Function g should satisfy the following equality for fractional area (an ) and volume (vn ):

∫ HR
HL
HR

∫

gdh = an,

(17.80)hgdh = vn .

(17.81)HL

We adopt a linear function of thickness for g. The thickness space is transformed to η = h − H L , and the thickness
distribution function is written as g = g1 η + g0 . These are substituted into (17.80) and (17.81) to yieldη 2R+ g0 η R = an,2
η3
η2
g1 R + g0 R = an η n,
3
2
g1

(17.82)(17.83)where η R = HR − H L and η n = hn − H L . These are algebraically solved for g0 and g1 as)
6an ( 2η R
−
η
,
n
3
η 2R
12an (
ηR )g1 = 3 η n −
.
2
ηR

g0 =

(17.84)(17.85)The values of the thickness distribution function at category boundaries are given as follows:)
6an ( 2η R
−
η
,
n
3
η 2R
6an (
ηR )
.
g(η R ) = 2 η n −3
ηRg(0) =(17.86)(17.87)Equation (17.86) gives g(0) < 0 when the thickness is in the right third of the thickness range or η n > 2η R /3. Equation(17.87) gives g(η R ) < 0 when the thickness is in the left third of the thickness range or η n < η R /3. Since a negative
g is physically impossible, we redefine the range of the thickness distribution function. Specifically, when the thickness
is within the left third of the thickness range, a new right boundary is set at HC = 3hn − 2H L and g is set to zero for[HC , HR ]. In this case, η R = HC − H L in (17.84) and (17.85). When the thickness is within the right third of the thickness
range, a new left boundary is set at HC = 3hn − 2HR and g is set to zero for [H L, HC ]. In this case, η R = HR − HC andη n = hn − HC in (17.84) and (17.85).
Finally, we remap ice into the original categories using the above thickness distribution function. If Hn∗ > Hn , ice istransferred from category n to n + 1. The transferred area ∆an and volume ∆vn are
∆an =

∫ Hn∗
gdh
Hn

– 206 –

(17.88)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea ice

and

∫ Hn∗

∆vn =

hgdh.

(17.89)Hn

If Hn∗ < Hn , ice is transferred from category n + 1 to n. The transferred area ∆an and volume ∆vn are
∆an =
and

∫ Hn
gdh

(17.90)hgdh.

(17.91)

Hn∗

∫ Hn

∆vn =

Hn∗

Snow and thermal energy are also transferred in proportion to the transferred volume. For example, ∆vsn = vsn (∆vin /vin )for snow and ∆ein = ein (∆vin /vin ) for thermal energy.
If ice is created in open water, the left boundary of category 1 (H0 ) is moved to F0 ∆t, where F0 is the growth rate inopen water. After area and volume are remapped in higher categories, ice area, volume, and energy are added to category1.
If ice is not created in open water, H0 remains zero, but the growth rate at the left boundary of category 1 is set toF0 = f 1 . If F0 < 0, the fractional area of category 1 thinner than ∆h0 = −F0 ∆t is added to open water area. In thisoperation, volume and energy are invariant. The area to be added to open water is
∆a0 =

∫ ∆h0
gdh.

(17.92)0

The right boundary of the thickest category N (H N ) is a function of its mean thickness h N . When h N is given, H N iscomputed as H N = 3h N − 2H N −1 . It is guaranteed that g(h) > 0 for H N −1 < h < H N and g(h) = 0 for H N < h.

17.4 Dynamics
17.4.1 Momentum equation for ice packThe momentum equation for an ice pack with mass ρ I Ah I is
1 ∂h∂
( Ah I u I ) − ρ I Ah I f vI = − ρ I Ah I g+ Fµ (σ) + A(τAI x + τIO x ),∂thµ ∂ µ1 ∂h∂ρ I ( Ah I vI ) + ρ I Ah I f u I = − ρ I Ah I g+ Fψ (σ) + A(τAI y + τIOy ),
∂t
hψ ∂ψ

ρI

(17.93)(17.94)

where (u I , vI ) is the velocity vector, h is the sea surface height, (Fµ, Fψ ) is the ice’s internal stress (which is a function of−−→internal stress tensor (σ)), and −τ→AI and τIO are stresses exerted by the atmosphere and ocean.

17.4.2 Stresses at top and bottom
The stress at the top is wind stress:−
τ→AI = CDAI ρ a |U a − u I |[(U a − u I ) cos θ a + k × (U a − u I ) sin θ a ],(17.95)where Ua is the surface wind vector, CDAI is the bulk transfer coeﬃcient between air and ice, ρ a is the density of air, andθ a is the angle between the wind vector and the ice drift vector, which is set to zero in MRI.COM.
Stress at the bottom is ocean stress:−τ−→
IO = CDIO ρo |Uw − u I |[(Uw − u I ) cos θ o + k × (Uw − u I ) sin θ o ],(17.96)

where Uw is the velocity of the first level of the ocean model, CDIO is the bulk transfer coeﬃcient between the ice andocean, ρo is the density of sea water, and θ o is the angle between the ice drift vector and the surface velocity of the ocean,which is set to 25° in the northern and −25° in the southern hemisphere following Hunke and Dukowicz (2002).

– 207 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.4 Dynamics

17.4.3 Internal stressIn a highly concentrated icepack, the eﬀect of the internal stress is as large as the Coriolis eﬀect and the surface stresses.
The expression of the internal stress is derived by regarding the ice as continuous media. The elastic-plastic-viscous (EVP)
model by Hunke and Dukowicz (1997, 2002) is adopted for the constitutive law (the relation between stress and strainrate). The EVP model is a computationally eﬃcient modification of the viscous-plastic (VP) model (Hibler, 1979). In the
VP model, the internal stress could be very large when the concentration is high and strain rate is near zero, which makes
the explicit integration infeasible. An alternative, the implicit method, is usually adopted, but it is not suitable for parallel
computing. The EVP model treats the ice as an elastic medium and a large local force is released by elastic waves, which
would be damped within the time scale of the wind forcing.
The constitutive law of the EVP model is1η−ζP
1 ∂σi j+
σi j +
σkk δi j +δi j = ϵ˙i j , i, j = 1, 2,
E ∂t
2η
4ηζ
4ζ

(17.97)

where ζ and η are viscous parameters, P represents ice strength, and E is an elastic parameter (mimics Young’s modulus).
In the VP model, tendency terms are zero.
The r.h.s. (ϵ˙i j ) is the strain rate tensor, expressed in Cartesian coordinates as:
ϵ˙i j =

1 ( ∂u I i ∂u I j )
+
.
2 ∂xj
∂ xi

(17.98)The divergence, tension, and shear of the strain rate are defined as follows:
DD = ϵ˙11 + ϵ˙22, DT = ϵ˙11 − ϵ˙22, DS = 2ϵ˙12 .

(17.99)The equation for the stress tensor for σ1 = σ11 + σ22 and σ2 = σ11 − σ22 is given by1 ∂σ1 σ1
P
+
+
= DD,
E ∂t
2ζ 2ζ
1 ∂σ2 σ2+= DT ,E ∂t2η
1 ∂σ12 σ12 1+= DS .
E ∂t2η
2(17.100)(17.101)(17.102)

In generalized orthogonal coordinates, divergence, tension, and shear of the strain rate are expressed by
1 [ ∂(hψ u I ) ∂(hµ vI ) ]+,h µ hψ∂µ∂ψhψ ∂ ( u I ) h µ ∂ ( v I )DT =−,h µ ∂ µ hψhψ ∂ψ hµh µ ∂ ( u I ) hψ ∂ ( v I )DS =+.
hψ ∂ψ hµh µ ∂ µ hψDD =(17.103)(17.104)(17.105)

The internal stress is obtained as the divergence of the internal stress tensor,)]
1 [ 1 ∂σ1
1 ∂(hψ σ2 )2 ∂ ( 2Fµ =++ 2
hµ σ12 ,2
2 hµ ∂ µ∂µ∂ψh µ hψ
h µ hψ

(17.106))]
1 [ 1 ∂σ1
1 ∂(hµ σ2 )2 ∂ ( 2− 2+
hψ σ12 .
2
2 hψ ∂ψhµ hψ ∂ψh µ hψ ∂ µ(17.107)

2

2

Fψ =

The viscous parameters are obtained from the ice strength and velocity as follows:P,2∆P
η= 2 ,
2e ∆
[
] 1/2
1
2
∆ = DD+ 2 (DT2 + DS2 ).
eζ=– 208 –

(17.108)(17.109)(17.110)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea ice

The pressure (strength) of the ice is a function of ice concentration and thickness:
P = P∗ Ah I e exp[−c∗ (1 − A)],

(17.111)

where P∗ is the scaling factor for pressure, c∗ is a parameter that defines the dependency on concentration, and e is theaxis ratio of the elliptic yield curve (e = 2).
The elastic parameter E is given by
2Eo ρ I Ah I
E=min(∆x 2, ∆y 2 ),(17.112)∆t e2where Eo is a tuning factor that satisfies 0 < Eo < 1, ∆t e is the time step for ice dynamics, and ∆x and ∆y are the zonal
and meridional grid widths.

17.4.4 Boundary conditions
Surface stresses on the ice are exerted for the fractional area of the ice within a grid cell. The ice concentration is multipliedby the wind and ocean stresses.
For the stress on the ocean, the ice-ocean stress is exerted for the ice-covered area, and the wind stress is exerted for the
open water area:( ∂U ∂V )A(1 − A)νV,= − (τIO x , τIOy ) +(τAO x , τAOy ).
(17.113)
k=0
∂z ∂z
ρo
ρo
Note that (τIO x , τIOy ) is reversed in sign because it is defined by (17.96) as the stress on the ice.

17.4.5 Solution procedure
Given the surface wind vector and the surface velocity of the ocean needed to compute surface stresses, the momentum
equations ((17.93) and (17.94)) and the equations for stress tensors ((17.100), (17.101), and (17.102)) are solved.
First, the stress tensor is computed using the equations for stress tensors, the momentum equation is then solved using
the stress tensor. Basically, the implicit method is used for prognostic variables for each equation. For example, stresstensor σ1 is solved for σ1m+1 as follows:
m+1
m
σ m+1
P
1 σ1 − σ1m+ 1 m + m = DD
.
E
∆t
2ζ
2ζ

(17.114)

Note that strain rate tensors and viscous parameters are updated every time step using a new velocity.
The momentum equations are solved using σ m+1 above:
ρ I Ah I

um+1
− um
1 ∂h
I
I
= ρ I Ah I f vIm+1 − ρ I Ah I g+ Fµ (σ m+1 ) + AτAI x∆thµ ∂ µ(17.115)

m+1+ Acw ρo |Uw − um
) cos θ o + (Vw − vIm+1 ) sin θ o ],I |[(Uw − u I

ρ I Ah I

vIm+1 − vIm
1 ∂h= − ρ I Ah I f um+1
− ρ I Ah I g+ Fψ (σ m+1 ) + AτAI y
I
∆t
hψ ∂ψ

(17.116)

m+1+ Acw ρo |Uw − um
) cos θ o − (Uw − um+1
) sin θ o ].
I |[(Vw − v I
I

Note that the surface velocity of the ocean (Uw, Vw ) is constant during the integration. The starting time level of the ocean
model is used, n − 1 for the leap-frog time step, and n for the Matsuno scheme. For the leap-frog time step of the ocean
model,
−
n−1
m
n−1
m+1
n−1
τ−→
) cos θ o + k × (Uw− um+1) sin θ o ].
(17.117)
IO = cw ρo |Uw − u I |[(Uw − u I
I
The time step of the ice dynamics is limited by the phase speed of the elastic wave. To damp the elastic waves during
the sub-cycle, the ice dynamics is sub-cycled several tens of steps during one ocean model time step.

– 209 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.5 Advection and Diﬀusion

17.5 Advection and Diﬀusion
Fractional area, snow volume, ice volume, ice energy, and ice surface temperature (optional; set flg_advec_tskin =.true. in mod_seaice_cat.F90) of each category are advected. A multidimensional positive definite advection transport
algorithm (MPDATA; Smolarkiewicz, 1984) is used. For all quantities which are advected, the harmonic-type diﬀusion
may be applied to remove noises.
The advection-diﬀusion equation for a property (α) is given by}{∂(hψ uα) ∂(hµ vα)1∂α++= D(α),(17.118)∂th µ hψ
∂µ∂ψ
where D represents diﬀusion operator. The specific representation for α is A for fractional area, Ahs for snow volume,
Ah I for ice volume, and Ah I E for ice energy.
In MPDATA, (17.118) is first solved to obtain a temporary value (α ∗ ) using the upstream scheme with a mid-point
velocity between time levels n and n + 1. Using this temporary value, an anti-diﬀusive velocity is computed as[]) ∂α ∗ 1∗111 ∂α
1 1 ( n+ 1
ũi+ 1 , j = ∗|u 2 |∆x − ∆t(un+ 2 ) 2− ∆tun+ 2 v n+ 2
2
α 2
∂x
2
∂y
1

1

1
1
∂u n+ 2 ∂v n+ 2 +− ∆tun+ 2 *
+
,
2
∂x
∂y
,
[
]
)
∗
∗
1 1 ( n+ 1
1
n+ 21 2 ∂αn+ 21 n+ 12 ∂α2|∆y − ∆t(v)
ṽi, j+ 1 = ∗|v− ∆tu v
2
α 2
∂y
2
∂x
1

(17.119)

1

1
1∂u n+ 2 ∂v n+ 2 ++− ∆tv n+ 2 *
.
2
∂y
, ∂x

(17.120)

This velocity is used to compute a new value using the upstream scheme starting from the above temporary value. Since
MPDATA is positive definite, the new area and thickness should be positive. If the sum of the fractional area exceeds
one, the ridging scheme will adjust the fractional area. Since energy is negative definite, the sign is reversed just before
advection and returned to a negative value after the advection.
The specific form of the harmonic-type diﬀusion is represented as follows:
{
(
)
(
)}
1
∂ hψ κ H ∂α∂ hµ κ H ∂αD(α) =+(17.121)h µ hψ ∂ µhµ ∂ µ∂ψhψ ∂ψ∂(hψ Fµα ) ∂(hµ Fψα ) 


1 


=−
+
(17.122)h µ hψ  ∂ µ∂ψ where κ H is the horizontal diﬀusion coeﬃcients. Using the flux form (17.122), the diﬀusion term is advanced in time withthe advection term (at the "correction" step of the MPDATA scheme).

17.6 Ridging
As a result of advection, the sum of the fractional area might exceed one, especially where the velocity field is convergent.
In such a case, it is assumed that ridging occurs among ice to yield a sum equal to or less than one. Even if the sum is lessthan one, ridging or rafting might occur where the concentration is high.
The ridging scheme of MRI.COM follows that of CICE, which is briefly summarized in this section.
First, we determine a fractional area that undergoes ridging: a P (h) = b(h)g(h). The weighting function b(h) is chosensuch that the ridging occurs for thin ice:
G(h) )2 (1−
if G(h) < G∗
∗
G
G∗
=0
otherwise,b(h) =(17.123)(17.124)

where G(h) is the area of ice thinner than h and G∗ is an empirical constant with G∗ = 0.15. The participation function
for category n (a Pn ) is obtained by integrating a P (h) for a range [Hn−1, Hn ] asa Pn =(
2
G n−1 + G n )
(G
−
G
)
1
−
,
n
n−1
G∗
2G∗
– 210 –

(17.125)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea icewhere dG = gdh is used. The property a Pn is the fractional contribution from the category n among the total area of ice
subject to ridging. The property G n is the area summed from category 0 to n. This equation is used for the category thatsatisfies G n < G∗ . If G n−1 < G∗ < G n , then G∗ is replaced by G n . If G n−1 > G∗ , then a Pn = 0. If the fractional area ofopen water exceeds G∗ (a0 > G∗ ), then ridging does not occur.
Ridging occurs such that the total area is reduced while conserving ice volume√ and energy. It is assumed that ice ofthickness hn is homogeneously distributed between Hmin = 2hn and Hmax = 2 H ∗ hn after ridging, where H ∗ = 25 m.
The thickness ratio before and after the ridging is k n = (Hmin + Hmax )/(2hn ). Therefore, when an area of category n is
reduced by ridging at a rate r n , the area of thicker categories is increased by r n /k n .
Among the new ridges, the fractional area that is distributed in category m is:
area
fm
=

HR − HL,
Hmax − Hmin

(17.126)

where H L = max(Hm−1, Hmin ) and HR = min(Hm, Hmax ). The fractional volume that is distributed in category m is:
vol
fm
=

(HR ) 2 − (H L ) 2.
(Hmax ) 2 − (Hmin ) 2(17.127)

The snow volume and ice energy are distributed by the same ratio as the ice volume.
The net area lost by ridging and open water closing is assumed to be a function of the strain rates. The net rate of area
loss of the ice pack (Rnet ) is given by
Rnet =

Cs(∆ − |DD |) − min(DD, 0),2(17.128)

where Cs is the fraction of shear dissipation energy that contributes to ridge building (0.5 is used in MRI.COM), DD is][
2 + 1 (D 2 + D 2 ) 1/2 . These strain rates are computed by the dynamics scheme.
the divergence, and ∆ = DD
2
T
S
e
∑N
The total rate of area loss due to ridging, Rtot = n=0
r n , is related to the net rate as follows:[

Rnet = a P0 +

N
∑

(
1 )]
Rtot .
a Pn 1 −kn
n=1(17.129)

Since Rnet is computed from (17.128), Rtot is computed from (17.129). Thus, the area subjected to ridging from categoryn is computed as ar n = r n ∆t (r n = a Pn Rtot ). The area after ridging is ar n /k n and the volume after ridging is ar n hn .
Using these, the ice subjected to ridging is first removed from category n, and the ridged ice is then redistributed into eachcategory.
In practice, we require that ar n ≤ an . If A > 1 after ridging, Rnet is adjusted to yield A = 1, and the ridging procedure
is repeated.

17.7 Adjustment
As the final procedure in each time step, the following adjustment is applied to the sea ice distribution.
• If the thickness of sea ice of some category is out of its category bounds, the sea ice in that category is moved tothe appropriate category.
• If the thickness of sea ice in category 1 is less than 0.1 meter, the thickness is set to 0.1 meter and the area fractionis adjusted accordingly.
• If the fractional area that sea ice occupies is too small, the sea ice of that category is forcibly melted (See Table
17.8).

17.8 Discretization
17.8.1 Advection (MPDATA)In MPDATA, tracer (α in Figure 17.3) is updated following a three-step procedure.

– 211 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.8 Discretization

a. A temporary value is computed using an upstream scheme.

The tracer fluxes at the side boundaries of a T-cell are:
] ∆y
(
) [
n
n
n
n
n
n
n
n
Fx αi,n j , αi+1,
,
u
+
|u
|)
+
α
(u
−
|u
|)
=
α
(u
,
i, j i+ 1 , j
j i+ 1 , j
i+1, j i+ 1 , j
i+ 21 , ji+ 21 , j
2
2
2
2
[
]
)
(
∆x
,
Fy αi,n j , αi,n j+1, vi,n j+ 1 = αi,n j (vi,n j+ 1 + |vi,n j+ 1 |) + αi,n j+1 (vi,n j+ 1 − |vi,n j+ 1 |)
2
2
2
2
2
2

(17.130)

where
)
1(
ui+ 1 , j+ 1 + ui+ 1 , j− 1 ,
2
2
2
2
2
2
)
1(
vi, j+ 1 = vi+ 1 , j+ 1 + vi− 1 , j+ 1 .
2
2
2
2
2
2

ui+ 1 , j =(17.131)(17.132)

The zonal flux is defined at the closed circle and the meridional flux is defined at the closed square in Figure 17.3.
Using this, the temporary value (α ∗ ) is computed using an upstream scheme:(αi,∗ j − αi,n j )∆Si, j∆t((n+ 21 )
n+ 21 )nnn
− Fx αi,n j , αi+1,= Fx αi−1,j, u 1
j , α i, j , u 1
i− 2 , j(i+ 2 , j(n+ 1 )
n+ 1 )+ Fy αi,n j−1, αi,n j , v 2 1 − Fy αi,n j , αi,n j+1, v 2 1 .
i, j−
i, j+
2

(17.133)

2

b. Compute the anti-diﬀusive transport velocity.

Using the temporary value computed in the first step, the anti-diﬀusive transport velocity is computed as follows:


)
∗
∗

1  1 ( n+ 12
1
n+ 21 2 ( ∂α )n+ 12 ( n+ 1 ∂α )ũi+ 1 , j =  ∗(2) |u 1 |∆x − ∆t(u 1 )− ∗(6) ∆tu 1 v 2
1
1 
i+
i+
,
j
,
j
,
j
,
j
,
j
i+
i+
i+
2
2 α 1
∂x
∂y
2
2 
2
2
2
α 1
 i+ 2 , j
i+ 2 , j1()∂u ∂v n+ 2
1
n+ 1,
− ∆tu 12+i+ 2 , j ∂ x
2∂ y i+ 12 , j


)
∗
∗

1
1  1 ( n+ 21
n+ 21 2 ( ∂α )n+ 12 ( n+ 1 ∂α )ṽi, j+ 1 =  ∗(2) |v 1 |∆y − ∆t(v 1 )−∆tv 1 u 2i, j+ 2i, j+ 2i, j+ 222 α 1∂ y i, j+ 12 α ∗(6) 1∂ x i, j+ 12  i, j+ 2
i, j+ 2() n+ 12
1
∂u ∂v
n+ 1
− ∆tv 2 1,+
i, j+ 2 ∂ x
2∂ y i, j+ 1

(17.134)(17.135)

2

where
1
∗
α ∗(2)1 = (αi,∗ j + αi+1,j ),i+ 2 , j21
α ∗(2) 1 = (αi,∗ j + αi,∗ j+1 ),i, j+ 2
2
1
∗
∗
∗
∗
∗
α ∗(6)1 = (αi,∗ j−1 + αi+1,
j−1 + 2α i, j + 2α i+1, j + α i, j+1 + α i+1, j+1 ),i+ 2 , j
8
1 ∗
∗
∗
∗
∗
∗
α ∗(6) 1 = (αi−1,
j + α i−1, j+1 + 2α i, j + 2α i, j+1 + α i+1, j + α i+1, j+1 ),i, j+ 28( ∂α ∗ )

=

∗
∗
αi+1,
j − α i, j,
∂ x i+ 21 , j∆xαi,∗ j+1 − αi,∗ j( ∂α ∗ )=,∂ y i, j+ 12∆y– 212 –

(17.136)(17.137)(17.138)(17.139)(17.140)(17.141)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea ice∗∂α ∗ )
1 [ n+ 1 ( ∂α ∗ )n+ 12 ( ∂α )= v 211
1 + v
1i, j− 2 ∂ y i, j− 21∂ y i+ 2 , j 4 i, j+ 2 ∂ y i, j+ 2∗∗]
n+ 1 ( ∂α )n+ 12 ( ∂α )+ v 211 + v
1 ,1i, j+ 2 ∂ y i+1, j+ 2
i, j− 2 ∂ y i+1, j− 2∗
[ n+ 1 ( ∂α ∗ )(1 ( ∂α ∗ )1 ∂α )
1
n+
un+ 2
= u 2
+ u 12i− 2 , j ∂ x i− 12 , j∂ x i, j+ 12 4 i+ 12 , j ∂ x i+ 21 , j( ∂α ∗ )( ∂α ∗ )
]
n+ 1
n+ 1
+ u 12+ u 12,1
1
i+ 2 , j+1 ∂ x i+ 2 , j+1
i− 2 , j+1 ∂ x i− 2 , j+1(

1

v n+ 2[]∂v ) n+ 21
1 ( ∂u ∂v ) n+ 12 ( ∂u ∂v ) n+ 12
=
+
+
+
,
∂ x ∂ y i+ 12 , j 2 ∂ x ∂ y i, j
∂ x ∂ y i+1, j[]( ∂u ∂v ) n+ 21
1 ( ∂u ∂v ) n+ 12 ( ∂u ∂v ) n+ 12
,
+
=
+
+
+
∂ x ∂ y i, j+ 12 2 ∂ x ∂ y i, j
∂ x ∂ y i, j+1
( ∂u

and

(

+

(17.142)(17.143)(17.144)(17.145)( n))n
u 1 ∆y i+ 1 , j − un 1 ∆y i− 1 , j + v n 1 ∆x i, j+ 1 − v n 1 ∆x i, j− 1
i+ 2 , ji− 2 , ji, j+ 2i, j− 2
∂u ∂v
2
2
2
2
+
=
.
∂ x ∂ y i, j
∆Si, j(17.146)c. Update tracer using the anti-diﬀusive velocity starting from the temporary value.
∗(αi,n+1
j − α i, j )∆Si, j

∆t

)
(
)
(
∗
∗
∗
∗
= Fx αi−1,
j , α i, j , ũi− 12 , j − Fx α i, j , α i+1, j , ũi+ 12 , j()()+ Fy αi,∗ j−1, αi,∗ j , ṽi, j− 1 − Fy αi,∗ j , αi,∗ j+1, ṽi, j+ 1 .
2

Ui-1/2, j+3/2

α i+1, j+1

α i, j+1

Ui+3/2, j+1/2

Ui+1/2, j+1/2
α i, jUi-1/2, j-1/2(17.147)

Ui+3/2, j+3/2

Ui+1/2, j+3/2

Ui-1/2, j+1/2

2

α i+1, j
Ui+1/2, j-1/2

Ui+3/2, j-1/2

Figure17.3 Position for tracer (α) and velocity (U). Area and thickness are defined at tracer points. Zonal fluxes arecomputed at closed circles and meridional fluxes are computed at closed squares. The budget is computed for a unit
cell for α.

– 213 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.9 Technical Notes

17.8.2 Momentum equation
Specific forms of discretization for properties related to internal stress are given here.
The strain rate tensor (ϵ)˙ and stress tensor (σ) are defined at tracer points (Figure 17.4).
Components (divergence, tension, and shear) of the strain rate tensor are expressed in a discretized form as follows:∆y i− 1 , j( ∆y i+ 21 , j
1
2(u I i+ 1 , j+ 1 + u I i+ 1 , j− 1 ) −(u I i− 1 , j+ 1 + u I i− 1 , j− 1 )
2
2
2
2
2
2
2
2
∆x i, j ∆yi, j22∆x i, j+ 1
∆x i, j− 1)
2
2+
(vI i+ 1 , j+ 1 + vI i− 1 , j+ 1 ) −(vI i+ 1 , j− 1 + vI i− 1 , j− 1 ) ,
2
2
2
2
2
2
2
2
2
2

(DD )i, j =(DT )i, j =u I i− 1 , j+ 1 ) ∆yi, j− 1 ( u I i+ 1 , j− 1
u I i− 1 , j− 1 )]1 [ ∆yi, j+ 21 ( u I i+ 12 , j+ 12
2
2
2
2
2
2
2
−
+
−
2 ∆x i, j+ 1 ∆y i+ 1 , j+ 1 ∆y i− 1 , j+ 1
∆x i, j− 1 ∆y i+ 1 , j− 1 ∆y i− 1 , j− 1
2

−

2

2

2 ∆yi+ 1 , j ∆x i+ 1 , j+ 1
2(DS )i, j =2

1 [ ∆x i+ 21 , j ( vI i+ 12 , j+ 12
2

−

2

2

2

vI i+ 1 , j− 1 )2

2∆x i+ 1 , j− 1
2

+

2

2

2

2

∆x i, j− 1 ( vI i− 1 , j+ 12

2

2

∆yi, j− 1 ∆x i− 1 , j+ 1
2

2

−

2

2

vI i− 1 , j− 1 )]2

2

∆x i− 1 , j− 12,

2

vI i− 1 , j+ 1 ) ∆yi, j− 1 ( vI i+ 1 , j− 1
vI i− 1 , j− 1 )]1 [ ∆yi, j+ 21 ( vI i+ 12 , j+ 12
2
2
2
2
2
2
2
−
+
−
2 ∆x i, j+ 1 ∆y i+ 1 , j+ 1 ∆y i− 1 , j+ 1
∆x i, j− 1 ∆y i+ 1 , j− 1 ∆y i− 1 , j− 1
2

+

2

2

2

1 [ ∆x i+ 21 , j ( u I i+ 21 , j+ 12
2 ∆yi+ 1 , j ∆x i+ 1 , j+ 1
2

2

2

−

2

2

u I i+ 1 , j− 1 )2

2∆x i+ 1 , j− 1
2

+

2

2

2

2

∆x i, j− 1 ( u I i− 1 , j+ 12

2

2

∆yi, j− 1 ∆x i− 1 , j+ 1
2

2

2

−

2

u I i− 1 , j− 1 )]2

2

∆x i− 1 , j− 1
2.

2

The internal stress is defined at velocity points and computed from stress tensor as follows:(Fµ )i+ 1 , j+ 1 =2

2

1 [ 1 ( (σ1 )i+1, j+1 + (σ1 )i+1, j − (σ1 )i, j+1 − (σ1 )i, j )2 2∆x i+ 1 , j+ 1
2

1(
+
2

∆y 2i+1, j+ 1 [(σ2 )i+1, j+1 + (σ2 )i+1, j ] − ∆y 2i, j+ 1 [(σ2 )i, j+1 + (σ2 )i, j ] )22∆y 2i+ 1 , j+ 1 ∆x i+ 1 , j+ 1
2

+

2

2

2

2

2

2
2
( ∆x i+ 12 , j+1 [(σ12 )i+1, j+1 + (σ12 )i, j+1 ] − ∆x i+ 12 , j [(σ12 )i+1, j + (σ12 )i, j ] )∆x 2i+ 1 , j+ 1 ∆y i+ 1 , j+ 1
2(Fψ )i+ 1 , j+ 1 =(17.148)

2

2

2

,

2

1 [ 1 ( (σ1 )i+1, j+1 + (σ1 )i, j+1 − (σ1 )i+1, j − (σ1 )i, j )2 2∆y i+ 1 , j+ 1
2

2

∆x
[(σ )
+ (σ2 )i, j+1 ] − ∆x i+ 1 , j [(σ2 )i+1, j + (σ2 )i, j ] )1 ( i+ 12 , j+1 2 i+1, j+1
2
2

−

2

2

∆x 2i+ 1 , j+ 1 ∆y i+ 1 , j+ 1
2+

(17.149)

2

2

2

2
2
( ∆y i+1, j+ 12 [(σ12 )i+1, j+1 + (σ12 )i+1, j ] − ∆y i, j+ 21 [(σ12 )i, j+1 + (σ12 )i, j ] )∆y 2i+ 1 , j+ 1 ∆x i+ 1 , j+ 1
2

17.9 Technical Notes
17.9.1 Source codesa. Program files

The sea ice part consists of the following programs.
– 214 –

2

2

2

.

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea iceUi-1/2, j+3/2

Ui+1/2, j+3/2

σi, j+1,εi, j+1
Ui-1/2, j+1/2

σi+1, j+1,εi+1, j+1
Ui+3/2, j+1/2

Ui+1/2, j+1/2
σi, j,εi, j

Ui-1/2, j-1/2

Ui+3/2, j+3/2

σi+1, j,εi+1, j
Ui+1/2, j-1/2

Ui+3/2, j-1/2

Figure17.4 Position of variables used by dynamics scheme.

ice_param.F90:

defines basic constants (This file is created by running the script configure on the basis of
the template ’ice_param.F90.in’. The number of thickness categories (NUM_ICECAT) mustbe specified in configure.in before running the script.)
ice_aggr_vars.F90: defines variables averaged over all categories for use in dynamics and ocean model
ice_cat_com.F90:
defines common variables and arrays
ice_cat_albedo.F90: computes surface albedo and fraction of shortwave penetrating into ice
ice_cat_bulk.F90:
computes surface fluxes using bulk formula
ice_cat_nest.F90:
contains subroutines relevant to nesting
ice_date.F90:
computes date and time
ice_dyn.F90:
computes ice dynamics
ice_flux.F90:
computes air-ice interface processes
ice_grid.F90:
sets grid cells (substitution from ocean model)
ice_hist.F90:
computes and outputs history
ice_main_cat.F90:
calls subroutines (main part of the ice model)
ice_mpdata.F90:
computes advection term using MPDATA
ice_remapv.F90:
remaps ice into thickness categories after thermodynamics
ice_restart.F90:
reads and writes restart
ice_ridge.F90:
computes ridging
ice_time.F90:
manages calendar
mod_seaice_cat.F90: computes thermodynamics and various adjustment processes

b. Structure

The ice model can be only used as a part of the ocean model. Among the subroutines included in the above programs, only
subroutines in ice_main_cat.F90 are called from surface flux routines of the ocean model. Relation with the OGCM is
as follows.
ogcm__ini (ogcm.F90)|+ surfflux_ctl__ini (surfflux_ctl.F90)|+ si_initialize (ice_main_cat.F90)... initializes the ice model by reading parameters and creating grid cells
ogcm__prepare (ogcm.F90)– 215 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.9 Technical Notes|+ surfflux_ctl__read_restart (surfflux_ctl.F90)|+ si_read_restart (ice_main_cat.F90)... reads restart and sets initial stateogcm__run (ogcm.F90)|+ LOOP_TIME_STEP|||+ part_1|||||+ surfflux_ctl__main (surfflux_ctl.F90)
|
|
|
|
|
+ surfflux_seaice (surfflux.F90)
|
|
|
|
|
+ simain (ice_main_cat.F90)||... the main program and calls subroutines of the ice model|||+ surfflux_ctl__write_restart (surfflux_ctl.F90)|||+ si_write_restart (ice_main_cat.F90)!
... writes restart|+ surfflux_ctl__run_end (surfflux_ctl.F90)|+ si_run_end (ice_main_cat.F90)... clears clock

Other things to be noted are:
• The sea ice model uses the forward scheme in time integration. The sea ice model is not called in the backward
part of the Matsuno (Euler backward) scheme.
• History of the basic state (integrated or averaged over all categories) might be monitored by the history routine of
the OGCM (specify namelist nmlhs_ice and nmlhs_iceuv).

17.9.2 Coupling with an atmospheric modela. General features

In coupled mode (SCUPCGCM), the boundary between the atmospheric component and the ocean-ice component is at the
air-ice(snow) boundary. The fluxes above the air-ice boundary are computed by the atmospheric component and passed
to the ocean-ice component via the coupler Scup (Yoshimura and Yukimoto, 2008). All the information needed by
the oceanic component is received by calling cgcm_scup__get_a2o (cgcm_scup.F90) from ogcm__run (ogcm.F90)
at the beginning of a coupling cycle. The information needed by the ice part is extracted by calling get_fluxi_a2o(get_fluxs.F90) from iaflux (ice_flux.F90). The main part of the ice is solved using surface fluxes and ice surface
temperature from the atmospheric component.
In the atmospheric component, temperature in the atmospheric boundary layer and at the ice surface (T3 (tsfcin)) arecomputed along with heat flux in snow layer (Q S = QI2 (fheatu)) using ice temperature (T1 (t1icen)), snow thickness(hs (hsnwn)), and ice thickness (h I (hicen)) given by the oceanic component (see Figure 17.1 and Table 17.1). The
properties needed by the atmospheric component are sent via cgcm_scup__put_o2a (cgcm_scup.F90) from ogcm__run(ogcm.F90) at the end of a coupling cycle.
To conserve heat and water in the coupled atmosphere-ocean system, globally integrated heat and fresh water flux mustmatch between the atmospheric and oceanic components. To achieve this, surface fluxes are adjusted in several steps,
which are explained in the following.

– 216 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea iceb. Correct errors due to interpolation
transFirst, to resolve errors that arise during the transformation process, a globally constant adjustment factor for heat flux qadj
is introduced so that
∫ (∑
ncat
)
trans
anini qn + qadj
dS = Qatm(17.150)
ocean,
n=0

where anini is the area fraction of sea ice for category n at the time of exchange, qn is the net heat flux for category n, andQatm
ocean is the net heat flux of the atmospheric model integrated over the ocean. The adjustment factor is given by
[
trans
qadj
= Qatmocean −

∫ ∑
ncat

]
anini qn dS /Socean,(17.151)n=0trans is added to the net longwave radiation flux.
where Socean is the global ocean area of the ocean model. qadj
trans
For fresh water, f adj is introduced so that

∫ (∑
ncat

evap

anini f n

)
trans
atm
+ f prcp + f rof + f adj
dS = Focean,(17.152)n=0
evap

where anini is the area fraction of sea ice for category n at the time of exchange, f n is the evaporation or sublimation for
atm is the net fresh water of the atmosphericcategory n, f prcp is the precipitation, and f rof is the continental run oﬀ, and Focean
model integrated over the ocean. The adjustment factor is given by
∫ (∑
ncat
) ]
[ atm
evap
trans
f adj = Focean −anini f n + f prcp + f rof dS /Socean,(17.153)n=0trans is added to either precipitation or evaporation according
where Socean is the global ocean area of the ocean model. f adj
to its sign.
This operation is done at subroutine adjust_fluxes_global of get_fluxs.F90.

c. Taking into account of the evolution of sea ice state

Second, as the ocean-sea ice system evolves during a coupling cycle, the area fraction of sea ice changes from the one at
the time of data exchange. Therefore, the global integral would be diﬀerent from the initial state. To resolve this diﬀerence,
evo is introduced so that
another adjustment factor qadj
∫ (∑
ncat

evo
anevo qn + qadj

)

dS =

∫ ∑
ncat

anini qn dS,(17.154)n=0

n=0

where anevo is the evolved area fraction of sea ice for category n. The adjustment factor is given as follows:
∫ ∑
ncat
ncat
[∫ ∑
]
evo
qadj
=
anini qn dS −anevo qn dS /Socean .
n=0(17.155)n=0trans + q evo is added to the longwave radiation flux.
In total, qadj
adj
LW
LW
trans
evo
qadj
= qatmn+ qadj+ qadjn(17.156)

evo is introduced so
For fresh water, because evaporation or sublimation may be aﬀected by the evolution of sea ice, f adj
that
∫ (∑
∫ ∑
ncat
ncat
)
evap
evap
evo
anevo f n + f adj
dS =anini f n dS,(17.157)n=0

n=0

where anevo is the evolved area fraction of sea ice for category n. The adjustment factor is given as follows:
∫ ∑
ncat
ncat
[∫ ∑
]
evap
evap
evo
anini f n dS −anevo f n dS /Socean .
f adj
=
n=0

n=0

evo is added to either evaporation or precipitation according to its sign.
f adjThis operation is done at subroutine adjust_fluxes_local of get_fluxs.F90.

– 217 –

(17.158)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.9 Technical Notes

d. Correct errors due to solution method

Third, discrepancies due to solution methods adopted at the air-ice interface (Section 17.2.3) are corrected. To ensure astable model integration, the ice surface temperature is set as an average of those of the atmospheric component and theocean-ice component. This new temperature T3o is used to replace the upward longwave flux of the adjusted heat flux
AI ) that contains the net longwave flux (q LW ). Because the longwave flux q LW contains theabove the air-ice interface (qadj
adj
adj
′

LW is replaced with q LW in the following manner,
upward flux due to T3a of the atmospheric model, qadj
adj
′

LW
LW
qadj
= qadj− ϵ σ(T3a ) 4 + ϵ σ(T3o ) 4 .

(17.159)′′

AI that contains q LW , theFor the surface where snow or ice is melting, the amount of melting is calculated using qadj
adj
discrepancy of the longwave radiation over melting surface should be stored for flux correction.
For the surface where melting does not occur, the interior flux (Q S or QI2 , they are equal), is obtained by equating′
AI
AI owing to the adjustment for the longwave radiation explained above. Furthermore,
qadj and Q S . Q S is diﬀerent than qadj′

AI and Q are not completely the same owing to the approximations used in the solution method explained in Section
qadj
S
17.2.3. This is not a serious problem for an ocean-only simulation. But this should be taken into account for the coupledmodel where exact surface heat flux conservation is required. Thus the diﬀerence between the adjusted air-ice flux and theAI − Q , should be stored.
interior flux, qadjS
To summarize, discrepancy of fluxes are stored in qdiﬀ in the following manner,
′

AI
AI
LW
LW
qdiﬀ = qadj
− qadj= qadj(T3a ) − qadj
(T3o ),
AI= qadj− Q S,

if surface is melting(17.160)

if surface is not melting.

(17.161)

stable , which
This is integrated over sea ice and divided by the oceanic area of the ocean model to give an oﬀset factor qadj
is added to the oceanic surface heat flux globally,

stable
qadj
=

∫ ∑
ncat

an (qdiﬀ )n dS/Socean .

(17.162)

n=1

Sampling is done at subroutine iaflux of ice_flux.F90 and flux adjustment is done at subroutine si_exit ofmod_seaice_cat.F90.

17.9.3 Nestinga. One-way nesting

In one-way nesting, it is recommended to set the boundary between coarse and fine resolution models (parent and childrespectively) along the line connecting tracer points. Setting velocity points as the boundary is awkward for several
reasons. For example, the velocity of the child model is zero at the boundary if sea ice of the parent model does not exist
around boundary. However, sea ice may exist inside the boundary in the child model in one-way nesting. Lacking thetransport removing sea ice from the child model through the boundary, sea ice may accumulate at tracer cells adjacent tothe boundary. A serious discontinuity may result. If tracer points are set as the boundary, this kind of accumulation will
not occur. However, conservation does not hold.
b. Two-way nesting

In two-way nesting, it is necessary to set velocity points as the boundary to impose global conservation including sea ice.
By reflecting the child model state to the parent model, the problem seen in one-way nesting, an inconsistent distribution ofsea ice around the boundary, may be avoided. However, owing to the application of the flux adjustment for conservation,
the positive definiteness of area fraction and volume of sea ice and snow, or negative definiteness of sea ice enthalpy are
not guaranteed after the advection equation is integrated.
This problem is avoided by the following treatment on the transport fluxes. First, the child model flux is adjusted tomatch the parent model flux (Chapter 18). Second, if the original flux calculated by the child model is going out of childmodel’s main region, the flux received from the parent model is replaced by this outgoing flux. Third, the diﬀerence
between the outgoing flux and the parent model flux is sent to the parent model. Fourth, if the value just inside theboundary of the child model is predicted to be negative, the boundary flux is adjusted so that the value does not get– 218 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea icenegative (a certain amount of quantity is taken from the parent). The adjustment factor is sent to the parent model. Fifth,in the parent model, the flux corrections received from the child model are summed over the corresponding parent gridand added to the original flux of the parent model.
Another problem arises for the interior temperature of sea ice. Owing to the discrepancy between the thermal energy
and the volume of the ice pack after advection, the interior temperature of sea ice layer may get erroneous. Therefore, onboth sides of the parent-child boundary, interior temperature of sea ice is calculated as the average of the top and bottom
surfaces of the sea ice,
T2 + T0adj
T1 =.
(17.163)
2
The diﬀerence between the predicted and adjusted ice energy is integrated along the parent-child boundary (Γ) as
Qdiﬀ =

∫ ∑
ncat
Γ n=1

[
]
adj adj
ρ I an hn E(T1 , r 1 ) − E(T1, r 1 ) dS.

(17.164)

This is reflected as a correction factor for the ocean surface heat flux of each model, which is computed asthermalcore
qadj= Qdiﬀ /Socean,(17.165)

core is the area of the main region of the parent or child model.
where Socean

17.10 Usage notes
17.10.1 Compilation
The information needed to compile the sea ice model with the OGCM should be given to configure.in. At least, ICE
option should be specified in the OPTIONS line and the number of categories of the sea ice thickness should be specified as
NUM_ICECAT, for example, a line NUM_ICECAT = 5 should be inserted somewhere in configure.in. The model options
related to the sea ice model are listed on Table 17.3.
Table17.3 Model options related to the sea ice model

option name
ICE
SIDYN

description
Sea ice model is included to the OGCM
Dynamics of sea ice is solved

CALALBSIHISTICECATAn albedo scheme of CICE is used
Monitor time evolution of basic state variables

ICEFULLMONIT

Activate extensive monitor

usage
otherwise, sea ice drifts with a third of the surface
ocean velocity.
otherwise constant
Use NAMELIST.name_model.MONITOR to specify theoutput group and interval. name_model is the nameof the model specified by configure.in
Use NAMELIST.name_model.MONITOR to specify theoutput group and interval

17.10.2 Job parameters (namelist)
The runtime job parameters (namelist) that should be given by NAMELIST.name_model are listed on Tables 17.4 to 17.11.
Job parameters related to monitoring are specified in the next subsection.
a. Thickness category

The bounds of thickness categories, whose number should be specified by NUM_ICECAT in configure.in before compilation, should be given at run time to namelist nml_seaice_thickness_category, which is explained in Table
17.4.

– 219 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.10 Usage notes
Table17.4 Namelist of thickness category of the sea ice model (nml_seaice_thickness_category)variable namehbound(0 : ncat)
lsicat_volchk

dimensionm
logicaldescription
category boundary
flag for checking mass conservation

usage (default value)
requireddefault = .false.

b. Time

The parameters related to time should be given by nml_seaice_time_conf and nml_seaice_time_calendar, which are
explained in Table 17.5. Note that the specification of nml_seaice_time_conf is not mandatory, necessary information
will be taken from the ocean model.
Table17.5 Namelist related to time of the sea ice model

variable namenstep_seaice
ibyri
ibmni
ibdyi
ibhri
ibmii
ibsci
l_force_leap

group
nml_seaice_time_conf
nml_seaice_time_conf
nml_seaice_time_conf
nml_seaice_time_conf
nml_seaice_time_conf
nml_seaice_time_conf
nml_seaice_time_conf
nml_seaice_calendar

l_leapnml_seaice_calendar

description
time steps to be proceedstart year of this run
start month of this run
start day of this run
start hour of this run
start minute of this run
start second of this run
forcibly specify leap year atrun time with l_leap
specification of whether thisyear is leap or not (valid whenl_force_leap = .true.)

usage (default value)same as OGCM
same as OGCM
same as OGCM
same as OGCM
same as OGCM
same as OGCM
same as OGCM
.false.
.false.

c. Dynamics

The parameters related to dynamics of sea ice (SIDYN option) should be given to namelist nml_seaice_dyn, which islisted on Table 17.6
Table17.6 Namelist related to dynamics of the sea icel (nml_seaice_dyn)

variable namedt_seaice_dyn_sec

dimension
secdescription
time step interval for dynamics (in seconds)

damp_vicem sec−1

damp_vmax

m sec−1

Rayleigh damping coeﬃcient for rapid velocity
minimum velocity for damping

usage (default value)
about a tenth of the baroclinic
time step of the ocean model
0.0
1.0

d. Diﬀusion

Coeﬃcient for the harmonic-type horizontal diﬀusion should be given to namelist nml_seaice_diff, which is explainedin Table 17.7.
Table17.7 Namelist of diﬀusion of the sea ice model (nml_seaice_diff)variable name
diff_seaice_m2ps

dimension
m2 sec−1

description
horizontal diﬀusion coeﬃcient

usage (default value)
requirede. Minimum fractional area

Minimum of fractional area applied to all thickness categories may be given at run time to namelist nml_seaice_model,
which is explained in Table 17.8.

– 220 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea ice
Table17.8 Namelist of minimum fractional area of the sea ice model (nml_seaice_model)

variable nameconcentration_min

dimension
1

description
minimum fractional area

usage (default value)
1 × 10−4

f. Air-ice transfer

The bulk transfer coeﬃcient at the air - sea ice interface should be given to nml_air_ice, which is listed on Table 17.9.
Table17.9 Namelist of air-ice bulk transfer coeﬃcient (nml_air_ice)variable name
chaiw
chait
chaie

dimension1
1
1

description
air-ice exchange coeﬃcient for windair-ice exchange coeﬃcient for sensible heat
air-ice exchange coeﬃcient for evaporation

usage (default value)
3.0 × 10−3
1.5 × 10−3
1.5 × 10−3

g. Restart

How to handle restart files is specified by namelists nml_seaice_run, nml_seaice_rst_in, and nml_seaice_rst_out,
which is explained in Table 17.10
Table17.10 Namelist related to restart of the sea ice model

variable namel_rst_seaice_in

group
nml_seaice_run

description

file_icecat_restart_in

nml_seaice_rst_in

nwrt_rst_ic

nml_seaice_rst_out

file_icecat_restart_out_temp

nml_seaice_rst_out

restart file name of categorized ice for input
the interval of time stepfor snap shot output
core part of the file namefor snap shot

usage (default value).true. : the initial stateis read from restart file.false.(default) : startfrom the state without sea
ice
required
required
requiredh. Albedo

When CALALBSI option is selected, an albedo scheme of CICE is used. The parameters of this scheme should be given bynamelist nml_albedo_seaice, which is listed on Table 17.11.
Table17.11 Namelist for the albedo of sea ice used by CALALBSI option (nml_albedo_seaice)variable namealb_ice_visible_t0
alb_ice_nearIR_t0
alb_snw_visible_t0
alb_snw_nearIR_t0
alb_ice_visible_dec_ratio
alb_ice_nearIR_dec_ratio
alb_snw_visible_dec_ratio
alb_snw_nearIR_dec_ratio

dimension(◦C) −1(◦C) −1(◦C) −1(◦C) −1

description
visible ice albedo for thicker ice
near infrared ice albedo for thicker ice
visible, cold snow albedo
near infrared, cold snow albedo
visible ice albedo declination rate
near infrared ice albedo declination rate
visible snow albedo declination rate
near infrared ice albedo declination rate

– 221 –

usage (default value)
0.78
0.36
0.98
0.70
0.075
0.075
0.10
0.15
Continued on next page

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.10 Usage notes

variable namehi_refatan_reftsfci_t0
tsfci_t1
fsnow_patch

Table 17.11 – continued from previous page
dimension
descriptionmthe maximum ice thickness to which connection function is used
the base value of the tangent hyperbolic connection function◦C
the temperature at which ice albedo isequated to that of ocean◦Cthe temperature at which ice albedo isstarted to decline to that of the ocean
meter
thickness of snow patch on melting bare ice

usage (default value)
0.50
4.0
0.0
−1.0
0.02

17.10.3 Monitoring
To monitor or sample the state of sea ice, HISTICECAT option must be specified. Extensive monitoring is activated byfurther specifying ICEFULLMONIT option.
The general behavior should be specified by the nml_seaice_budget and nml_seaice_hst in NAMELIST.name_model,
which is explained in Table 17.12.
Table17.12 Namelist of monitor of the sea ice model

variable namenstep_si_budget_intervalgroup
nml_seaice_budget

undef_value

nml_seaice_hst

undef_for_land

nml_seaice_hst

description
the interval of time step bywhich water budget is written
undefined value to fill ice freegrid point
undefined value to fill land
grid points

usage
0 : none, −1 : monthly
real(8)real(8)Sampling variables are grouped according to the processes they are involved, the grid points they are defined, and thepriorities given by Sea Ice Model Intercomparison Project (SIMIP) (Notz et al., 2016).
The groups are listed on Table 17.13. Users specify nml_history blocks in file NAMELIST.name_model.MONITOR in
the following manner.
&nml_history
namefile_baseinterval_step
suffix&end

Specification of NAMELIST.name_model.MONITORforseaicemonitoring= ’sea ice basic state vars at t-point’,= ’result/hs_ice_t_sta1’,= 12,= ’minute’,
Table17.13 List of groups for sea ice monitoringgroup namesea ice basic state vars at t-point
sea ice state vars at t-point
tendencies of sea-ice mass and area fraction at t-point
heat and freshwater fluxes over sea ice fraction at t-point
sea ice state vars of each category at t-point
sea ice basic dynamics at u-point

usage note"mean" is forced"mean" is forced"mean" is forced"mean" is forced"mean" is forced"mean" is forced
Continued on next page– 222 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 17

Sea ice

Table 17.13 – continued from previous pagegroup nameusage notesea ice dynamics at u-point"mean" is forced
sea ice dynamics at t-point"mean" is forced
sea ice vars at x-point
sea ice vars at y-point
sea ice integrated measures
sea ice dynamics at u-point priority 3
valid if ICEFULLMONIT, "mean" is forced
sea ice dynamics at t-point priority 3
valid if ICEFULLMONIT, "mean" is forced
sea ice fluxes of each category at t-point
valid if ICEFULLMONIT
heat and freshwater fluxes over open leads at t-point
valid if ICEFULLMONIT

17.11 Appendix
17.11.1 Partial derivative of specific humidity with respect to temperature
In MRI.COM surface temperature of ice is computed using semi-implicit method, where an expression for the partial
derivative of the saturation specific humidity with respect to temperature (t) is needed as in (17.35).
Properties of moist air used by MRI.COM is based on Gill (1982) and replicated in Section 14.12.3. Using (14.117),
the saturation specific humidity (qi ) is given asqi = ϵ esi′ /(ps − (1 − ϵ )esi′ ).

(17.166)

where saturation partial pressure of moist air over ice (esi′ ) is given by (14.123) and (14.124). Using (17.166),
∂esi′
∂qi
ϵ ps
=
,
′
2
∂t
{ps − (1 − ϵ )esi } ∂t(17.167)

where ∂esi′ /∂t is expressed by setting f w = 1 as∂esi′
= ln 10 · 10g(t) · g ′ (t),∂t(17.168)

where g(t) = (0.7859 + 0.03477t)/(1 + 0.00412t) + 0.00422t and g ′ (t) = ∂g(t)/∂t, where t is temperature at the uppersurface of sea ice in ◦C.

17.11.2 Physical constant, parameters
Since the ice part is coded in SI units, constants and parameters are written in SI units.
a. Thermodynamics
Table17.14 Physical constants and parameters relevant to sea ice thermodynamics

parameter
Thermal ice conductivity
Thermal snow conductivity
Specific heat of sea waterSpecific heat of air
Specific heat of icevaluek I = 2.04 J m−1 s−1 K−1
k s = 0.31 J m−1 s−1 K−1
Cpo = 3990 J kg−1 K−1Equation (14.125)
Cpi = 2093 J kg−1 K−1– 223 –

variable name in MRI.COM
cki
cks
cp0
cpair
cpi
Continued on next page

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201717.11 Appendix
Table 17.14 – continued from previous pagevalue
Specific heat of snow
Cps = 0.0 J kg−1 K−1
Stefan Boltzmann constantσ = 5.67 × 10−8 W m−2 K−4
Albedo of open ocean surface
α o = 0.1 (default)
Albedo of ice
α I = 0.6 (default)
Albedo of snow
α s = 0.75 (default)
Emissivity of ocean surface
ϵ o = 1.0
Emissivity of ice surface
ϵ I = 1.0
Emissivity of snow surface
ϵ s = 1.0
neutral bulk transfer coeﬃcient for heat
CHAI = 1.5 × 10−3 (default)
neutral bulk transfer coeﬃcient for moisture
CEAI = 1.5 × 10−3 (default)
Latent heat of fusion
L F = 3.347 × 105 J kg−1
Latent heat of sublimation
Equation (14.127)
Constants for fusion phasem = −0.0543 ◦C/pss
equation: T f = mS + nz
n = −0.000 759 ◦C m−1
Ice roughness parameterz0 I = 0.05h I /3
Salinity of sea ice
SI = 4.0 pss
von Karman’s constantk = 0.4
Thickness/compactness diﬀusion of iceκ H = 1.0 × 103 m2 s−1
Seawater kinematic viscosityν = 1.8 × 10−6 m2 s−1
Seawater heat diﬀusivity
α t = 1.39 × 10−7 m2 s−1
Seawater salinity diﬀusivity
α b = 6.8 × 10−10 m2 s−1
Turbulent Prandtl number
Prt = 0.85
b in eqs (17.57),(17.59)
b = 3.14
parameter

variable name in MRI.COM
—
stbl
albw
albi
albs
eew
eei
ees
chaih
chaiw
alf
rlth
xmxm
xnxn
z0
si
xk
akh
anu
at
as
prt
ab

b. Dynamics
Table17.15 Physical constants and parameters relevant to sea ice dynamics

parameter
Density of sea water
Density of air
Density of sea ice
Reference snow density (ratio between snow andwater)
e-folding constant for ice pressurepressure scaling factor
drag coeﬃcient (air-ice)drag coeﬃcient (ice-ocean)yield curve axis ratio
scaling factor for Young’s modulus
water turning angle
air turning anglevalueρo = 1036 kg m−3
ρ a , Equation (14.128)
ρ I = 900 kg m−3
ρs /ρo= 330 kg m−3 /1036 kg m−3c∗ = 20.0
P∗ = 2.75 × 104 N m−2
CDAI = 3.0 × 10−3 (default)
CDIO = 5.5 × 10−3
e = 2.0
Eo = 0.25
θ o = ±25° (positive/negative in thenorthern/southern hemisphere)θ a = 0° (not considered)– 224 –

variable name in MRI.COM
ro
roair
rice
rdsw
cstar
prsrefchaiw
cdrgiwelips
eyoungwiangl–TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 18

Nesting
In MRI.COM, embedding of a fine-resolution regional model within a coarse-resolution model can be realized by nesting
method. The available methods range from a serial execution of a coarse-resolution (parent) model writing a boundary data
and a fine-resolution (child) model reading it (oﬀ-line one-way nesting) to a parallel execution of both models exchanging
data in both directions (on-line two-way nesting), as summarized in Table 18.1.
This chapter explains facilities available in MRI.COM for nesting. How to construct and run a set of coarse- and fineresolution models in nesting are also explained. Usage of nesting in the sea ice model and the tidal scheme is described in
Section 17.9.3 and 6.5.3, respectively.

18.1 Feature
In a set of nested grid models, a fine-resolution (child) model is embedded in a coarse-resolution (parent) model. In
one-way nesting, values at the side boundary of the child model are given by the parent model. The side boundary data
may be transferred both oﬀ-line and on-line. In oﬀ-line mode, the data needed to calculate side boundary values are
written to files by first running the parent model, and the child model is executed reading these data and calculating the sideboundary values by interpolation. In on-line mode, parent and child models are run in parallel. Simple coupler (Scup),
which is originally developed by Yoshimura and Yukimoto (2008) for the communication among components of the MRI
Earth System Model, is used to transfer data.
In two-way nesting, in addition to the transfer of the side boundary data from parent to child model, the result of thechild model is reflected to the parent model in the embedded region. Data are exchanged on-line and Scup is used for datatransfer.
The boundary between parent and child models may be formed by connecting either the tracer points or the velocity
points of the Arakawa B-grid arrangement (Figure 18.1). Table 18.1 shows the side boundaries suited for the available
nesting methods. From experiences, we recommend setting the tracer points as the side boundary for one-way nesting(imposing Dirichlet boundary condition for the sea surface height and tracer equations) and the velocity points as the sideboundary for two-way nesting (imposing Neumann boundary condition for the sea surface height and tracer equations).
In one-way nesting, the child model is quite stable when the side boundary is set to the tracer points. In two-way nesting,
imposing the conservation on water volume and tracer content for the set of parent and child models is more straightforward
when the side boundary is set to the velocity (flux) points.
To be able to treat both side boundary conditions with a single topographic configuration, the domain of a child model
is constructed by taking tracer points as the boundary of the main region (Figure 18.1a). Inside the boundary, the same
topography as the parent model should be given for at least one parent grid cell. By doing so, we may use the velocity
points located inside the boundary by a half parent grid width as the parent-child boundary for two-way nesting (Figure
18.1b). The model domain needs additional two rows or columns of velocity cells outside the main region along each sideboundary (hereinafter referred to as halo; the region shaded by dark green in Figure 18.1a).
Whether the side boundary of the child model is the tracer points or the velocity points, three columns or rows of theparent model are suﬃcient for a child model to fill both halo and boundary regions. Here, the boundary region refers to a
Table18.1 Nesting methods available in MRI.COM and the recommended side boundary of a child model in eachmethod. It should be noted that the region and the topography of a child model must be made as if tracer-points areused as the side boundary.

model run
one-way nest
two-way nest

on-line (coupler)
tracer point
velocity point– 225 –

oﬀ-line (file I/O)
tracer point
unavailable

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201718.2 Grid and Topography
region inside the boundary of the main region by one grid of the parent model. For example, T1 − T3 and U1 − U3 in Figure18.1 are suﬃcient for filling the southern boundary (except for one-to-one nesting). In oﬀ-line one-way nesting, the parent
model outputs the three columns or rows to files, and the child model reads them. In on-line nesting, the mapping tablebetween parent and child grids should be prepared for the Scup library. The Scup library exchanges only necessary data
among MPI processes of diﬀerent models.
In this way, prognostic variables at the side boundary of the child model are directly replaced by those from the parent
model. This is usually called the clamped method (Cailleau et al., 2008). This method does not guarantee conservation of
tracers in contrast to the method where fluxes are imposed at the boundary, but the integration is quite stable.
(a)

(b)

U3

U3

T3
T2
T1

T3
U2

T2

U1

T1

×：tracer

U2
U1●：velocity

×：tracer●：velocity

Figure18.1 Relation between parent (large symbols) and child (small symbols) grid points. The region shaded withblue (green) is the main region of the child (parent) model. The boundary of the main region is indicated by the redlines. In (a), the tracer points, and in (b), the velocity points, are set as the boundary between the two models. The
domain of the child model should be constructed to include the dark green area, which is expanded outward by two
velocity cells from the tracer point boundary (red line in (a)). This expanded region is used as halos for the main region.
Variables in the halo region are obtained by using values of the parent model. Thus, for the southern boundary, threerows (T1 - T3 and U1 - U3 ) are sent to the child model from the parent model in both case (a) and (b). In two-way nesting,variables in the main region of the child model shown in (b) are sent to the parent model, except on the boundary line
(red line).

Multistage nesting is available in MRI.COM, and on-line nesting and oﬀ-line nesting can be used together. This
flexibility makes various kinds of experiments possible, e.g., a double on-line nesting experiment using global, NorthPacific, and Japanese coastal models, or an oﬀ-line one-way nesting experiment of a Japanese coastal model reading theside boundary values from an on-line nesting experiment for a set of global and North-Pacific models.

18.2 Grid and Topography
18.2.1 Parent (coarse-resolution) model
There is no particular issue in constructing the parent model. Because the grid size ratio between the child and parent
models must be odd, the grid size of the child model may be taken into account in determining the grid size of the parent
model. A grid size ratio of 1:3 or 1:5 is recommended.

18.2.2 Child (fine-resolution) modela. Horizontal grid size

The same grid size as the parent model should be used in the vertical direction. The child model is embedded in thehorizontal directions. The distance between the velocity points of the parent model (i.e., T-cell) is divided equally in thechild model.
– 226 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 18

Nesting

The domain of a child model is constructed by taking the tracer point as the boundary of the main region and adding
the halo region, even if the velocity-point boundary will be used. Here, the halo region is two rows or columns of velocity
cells outside the boundaries as shown in Figure 18.1a. There should be always two marginal velocity cells outside themain region at the western and southern boundaries, even if the western or southern boundary is filled by land and does
not receive data from the parent model. However, if the eastern or northern boundary is filled by land and does not receive
boundary data, one marginal velocity cell will be enough.
How to divide parent cells"dxt" of the parent model (U-point distance, T-point at its center) should be
divided equally in the child model.
There should be two marginal U-points in the child model(one cell will do if the eastern or northern boundary is filled by land)***** parent (coarse-resolution) model *****"ifst" (variable name is sub_region_ifirst) is the western end grid point numberof the parent model to be used for interpolation.
boundary|(i=ifst)(i=ifst)(i=ifst+1) (i=ifst+1)<------- dxt -------> U <------- dxt -------> U||-----------+-----------|-----------+-----------|---|
T <------- dxu -------> T <------- dxu -------> T(i=ifst)(i=ifst)(i=ifst+1) (i=ifst+1)|***** child (fine-resolution) model *****(nesting ratio is 1:3)|
(1) (1) (2) (2) (3) (3)<-dxt->U<-dxt->U<-dxt->U<-dxt->U|---+---|---+---|---+-T<-dxu->T<-dxu->T<-dxu->T(1) (1) (2) (2) (3) (3)b. Topography of the child model

To achieve flux conservation, the child model should have the same topography as the parent model around the sideboundary of the child model. It is recommended that the topography of the child model should have the same topography
as the parent model for two parent velocity cells inward and one parent velocity cell outward relative to the boundary of
the main region constructed by connecting tracer points (Figure 18.2).
c. Weighting ratio between parent and child models

For smoothness, prognostic values of the child model around the side boundary may be given as a weighted average of thesolutions of the parent and child models. For variable ϕ at a grid point with index (i, j),p

ϕchild= αi j ϕi j + (1 − αi j )ϕicj ,ij(18.1)

where ϕchild is the value for the child model, ϕ p is the parent model solution interpolated on the child grid, ϕc is the childmodel solution, and α is a spatially dependent weight between parent and child solutions.
Although this weight is introduced for a smooth transition from parent to child model solution, we recommend that theweight for the parent model to be unity in the halo and at the boundary of the main region and zero elsewhere. Instead,this weight is used to specify the boundary between the main regions of parent and child models. When the boundary isat the tracer points, the weight for the parent solution should be unity in the dark green region of Figure 18.1a. When theboundary is at the velocity points, the weight for the parent solution should be unity in the dark green region of Figure
18.1b. Weight must be prepared separately for the tracer and velocity points. (Incidentally, a radiation boundary condition
may be imposed by elaborating weight values.)– 227 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201718.3 Intermodel transfer

×：tracer●：velocity

Figure18.2 Condition on the topography of the child model around the side boundary. The topography of the childmodel should have the same topography as the parent model for two parent velocity cells inward (dark blue) and one
parent velocity cell outward (light blue) relative to the boundary of the main region constructed by connecting tracer
points (red rectangle).

18.3 Intermodel transfer
Communication between models is both downscaling and upscaling in two-way nesting, whereas only downscaling in
one-way nesting.

18.3.1 Downscaling
In the downscaling from parent to child model, the method of data transfer and interpolation is diﬀerent between tracer
point and velocity point.
a. Tracer points

It is straightforward to do a linear interpolation for variables on tracer points. Figure 18.3 shows the relative positionsof tracer points of parent and child models. Because topography is common between parent and child models aroundthe side boundary, all child tracer points in an ocean velocity cell of the parent model (bounded by the red rectangle) areocean points. Values on these child tracer points are computed by using values at the four corner tracer points of theocean velocity cell of the parent model (⊚ in Figure 18.3). All of the four corner points are always available for bilinear
interpolation to compute the tracer values of the child model. Thus, assuming that parent and child grid points coincide
at (i p, j p ) and (i c, jc ), the tracer value at (i c + m, jc + n), where (0 ≤ m ≤ M) and (0 ≤ n ≤ N ), and M and N are thezonal and meridional nesting ratios, is computed as follows:
p

p

p

p

ϕic +m, jc +n = (1 − w x )(1 − wy )ϕi p , j p + w x (1 − wy )ϕi p +1, j p + (1 − w x )wy ϕi p , j p +1 + w x wy ϕi p +1, j p +1,(18.2)

where w x = (x itc +m − x it p )/(x it p +1 − x it p ) and wy = (y tjc +n − y tj p )/(y tj p +1 − x tj p ), x t and y t represent the tracer point inzonal and meridional coordinates, respectively.
b. Velocity points

Figure 18.4 shows the relative positions of velocity points of parent and child models. Unlike the tracer points, all of thefour parent velocity points surrounding a child velocity point are not necessarily ocean points. When all parent grid points
are ocean, bi-linear interpolation can be applied.
A problem arises for velocity points of the child model between an ocean point and a land point of the parent model(shaded by blue in Figure 18.5a), because only data at a single point is available for interpolation. The same problem also
– 228 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 18

Nesting

Tic , jc+N

Tic+M, jc+N

Tic+M-1, jc+N

Tic+1, jc+N

Tip+1, jp+1

Tip, jp+1
Tic , jc+N-1

Tic+1, jc+N-1

Tic+M-1, jc+N-1

Tic+M, jc+N-1

Uip , jp
Tic -1, jc +1

Tic , jc +1

Tic+1, jc +1

Tic+M-1, jc+1

Tic+M, jc+1

Tic -1, jc

Tic, jc

Tic+1, jc

Tic+M-1, jc

Tic+M, jc

Tip , jp
Tic -1, jc -1

Tip+1, jp
Tic, jc -1

Tic+1, jc -1

Figure18.3 Relative positions of tracer points of parent and child models (⊚ and ◦, respectively), coinciding at (i p, j p )and (i c, jc ). The red and blue frames indicate U- and T-boxes in the parent model, respectively. M and N are the zonal
and meridional nesting ratios.

arises for a velocity cell that has a partial velocity cell at its neighbor. The velocity for the partially blocked part (shadedby purple in Figure 18.5a) should be treated separately from the upper unblocked part.
In this case, the child model velocity perpendicular to the coast line (uc ) is assumed to be zero at the coast and linearlyinterpolated from the nearest ocean grid point of the parent model (Figure 18.5b). The component tangent to the coast line
(v c ) is set to the same value at the nearest ocean grid point of the parent model (Figure 18.5c). By doing so, the interpolated
velocity field yields nearly the same transport along the coast line as the parent model. To be specific, when the parent
model has a partially shaved bottom cell at (i p − 1, j p, k), velocities of the child model between x u (i p − 1) ≤ x u ≤ x u (i p )are computed as follows:pp

uc (i) =w x ui p + (1 − w x )ui p −1,(18.3)ppv c (i) =w x vi p + (1 − w x )vi p −1,

uc (i) =v c (i) =∆z i p −1,k∆z i p ,k∆z i p −1,k∆z i p ,kp

if x u (i) < x t (i p )p∆z i p ,k − ∆z i p −1,kp∆z i p ,k∆z i p ,k − ∆z i p −1,k(w x ui p + (1 − w x )ui p −1 ) +p(w x vi p + (1 − w x )vi p −1 ) +∆z i p ,kp

wui p ,p

vi p ,(18.4)(18.5)

if x t (i p ) < x u (i)(18.6)

where w x = (x u (i) − x u (i p − 1))/(x u (i p ) − x u (i p − 1)) and w = (x u (i) − x t (i p ))/(x u (i p ) − x t (i p )).
To apply this method to the barotropic transport, the barotropic transport is first divided by the height of the watercolumn to obtain the vertically uniform barotropic velocity. This velocity is given to all the velocity points in the watercolumn and then the same method as used for the 3-D velocity field is applied.

– 229 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201718.3 Intermodel transfer

Uic , jc+N

Uic+1, jc+N

Uic+M-1, jc+N

Uip , jp+1

Uic -1, jc+1

Uip+1, jp+1

Uic , jc+N-1

Uic+1, jc+N-1

Uic+M-1, jc+N-1 Uic+M, jc+N-1

Uic , jc+1

Uic+1, jc+1

Uic+M-1, jc+1

Uic+1, jc

Uic+M-1, jc

Uic , jc

Uic -1, jc

Uic , jc -1

Uic+M, jc+1

Uic+M, jc
Uip+1, jp

Uip , jp
Uic -1, jc -1

Uic+M, jc+N

Uic+1, jc -1

Figure18.4 Same as Figure 18.3 except for velocity point.

18.3.2 Upscaling
In two-way nesting, values of the parent model in the child model’s main region are filled with those based on the childmodel. The region to receive feedback in the parent model is the main region of the child model excluding the boundaryline (inside the red line of Figure 18.1b). MRI.COM has two options for upscaling, replacement and nudging.
a. Replacement

One option is to simply replace variables of the parent model by those of the child model. Ideally, fine scale features thatcannot be resolved by the parent model should be removed in upscaling. In MRI.COM, simple averaging within a parent
grid is used (averaging over blue rectangle in Figures 18.3 and 18.4), though more sophisticated filters have been proposed(e.g., Debreu et al., 2012). Because applying a full averaging for the entire embedded region is costly in terms of datatransfer between models, a full averaging is only applied to the 2 to 3 grid points inside the parent-child boundary,p

ϕi p , j p =1 ∑ c
ϕ , i c − i m ≤ i ≤ i c + i m and jc − jn ≤ j ≤ jc + jn,
M × N i, j i, j

(18.7)

where i m = (M − 1)/2 and jn = (N − 1)/2. For the region further inside the boundary, values on the child grid points thatcoincide with the parent grid points may be used without filtering,p

ϕi p , j p = ϕicc , jc ,x t (i p ) = x t (i c ) and y t ( j p ) = y t ( jc ).

(18.8)

This treatment does not cause a problem because this inner region away from the boundary does not aﬀect the temporal
evolution of the physical field in the main region of the parent model.

– 230 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 18(a)Nesting
ip-1
ic -M

×：tracer / normal ﬂux

ip
ic

(b)

ip-1
ic -M

●：velocity

×：ﬂux

ip
ic

●：velocity

(c)

ip-1
ic -M

×：ﬂux

ip
ic

●：velocity

Figure18.5 Interpolation of the parent model velocity to the child grid with bottom topography. (a) For velocity pointsof the child model between an ocean point and a land point (shaded with blue or purple), only one data is available forinterpolation. Examples of (b) normal and (c) tangential component of velocity relative to the coast line (red arrows).
For an ocean velocity cell next to a partial cell, the cell is divided into a part blocked by a wall and a part unblocked(green point). Linear interpolation is used for an unblocked part and special treatment is applied to the blocked part
(purple arrows). The velocity is obtained as a weighted average of the two parts.

b. Nudging

Nudging is another option to reflect the child result to the parent model. In MRI.COM, tracer values of the child modelp
may be used as reference values (ϕC ) to which parent tracer fields (ϕ p ) are nudged,∂ϕ p1p= − (ϕ p − ϕC ),∂tγ(18.9)

where γ is a restoring time scale usually taken to be very small, about 0.1 day. Only tracer fields are corrected in the parent
model. Neither velocity nor sea surface height is modified. Our assumption lying behind this treatment is that reflecting
only slowly evolving quasi-geostrophic baroclinic fields is enough for a parent model in many applications.

18.4 Conservation
Conservation of sea water volume and scalar quantities are generally preferable. This is especially so for the case ofclimate applications with long-term integration. To achieve this in nesting, MRI.COM has options to adjust side boundaryand surface fluxes in the child model, or source / sink due to nudging in the parent model.

18.4.1 Side boundaries
Adjustment of the flux at the side boundary is performed as follows. Consider an evolution equation for a quantity ϕ in a
simple flux form:∂ϕ= −∇ · f,(18.10)∂twhere f = ( f x, f y ), which are due to either advective or diﬀusive fluxes. Integration (18.10) over each model’s main
region (Figure 18.6) gives
∫
∫
p
∂ΦΩ−ω
∂ϕ
≡
dS = −
f p · ndl
(18.11)
∂t
∂t
Ω−ω
∂Ω
for parent model’s main region and

c
∂Φω
≡
∂t

∫

∂ϕ
dS = −
ω ∂t– 231 –∫
∂ωf c · ndl

(18.12)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201718.4 Conservation
for child model’s main region. If integral values of the flux normal to the boundary (that is, r.h.s.) cancel between parent
and child models, the sum of integral values of the tendency of ϕ in each model’s main region is conserved,
p

∂ΦΩ−ω

+

∂t

c
∂Φω
= 0.
∂t

(18.13)

However, this conservation does not always hold owing to discretization. Some adjustment is needed.
In MRI.COM, conservation of water volume and tracers is achieved by adjusting boundary fluxes computed in eachmodel. Diﬀerent approaches are taken for advective and diﬀusive fluxes.
For advective fluxes, fluxes of the child model are adjusted to match those of the parent model. Fluxes for tracer cells of
the child model ( f yc ) are area-integrated along the side boundary section over the width of a tracer cell of the parent model(Figure 18.7a),∑f yc (i) × ∆s c (i), i c − i m ≤ i ≤ i c + i m,
(18.14)Fyc (i p ) =iwhere ∆s c is the area of the side boundary of a tracer cell of the child model and i m = (M − 1)/2. This is compared withp
the flux for the tracer cell of the parent model ( f y ) (Figure 18.7b),pp

Fy (i p ) = f y (i p ) × ∆s p (i p ),(18.15)where ∆s p is the area of the side boundary of a tracer cell of the parent model. On the basis of this comparison, fluxes of
the child model are adjusted using a quantity f yc adj (i p ) given byp

f yc adj (i p ) = (Fy (i p ) − Fyc (i p ))/∆S c (i p ),where
∆S c (i p ) =∑∆s c (i), i c − i m ≤ i ≤ i c + i m .

(18.16)(18.17)i

The adjusted flux f yc mod is given by
f yc mod (i) = f yc (i) + f yc adj (i p ), i c − i m ≤ i ≤ i c + i m,

(18.18)so that the flux integral in the child model equals the parent model flux for parent’s tracer cell.
For diﬀusive fluxes, the sum of fluxes in the child model ( f xc ) may be given to the parent model instead of adjusting thechild model flux to the parent model flux as in the advective fluxes (Figure 18.8),∑pp
Fx mod ( j p ) = f x mod ( j p )∆s p ( j p ) =f xc ( j) × ∆s c ( j), jc − jn ≤ j ≤ jc + jn .
(18.19)jp

p

Fx mod ( j p ) replaces the flux calculated on the parent grid Fx ( j p ) at the boundary.
The flux adjustment should be applied to any horizontal flux of volume and tracers. Currently, this is implemented in
the following fluxes:• Volume flux in the 3-D continuity equation• Advective volume flux in the 2-D continuity equation (both barotropic and baroclinic time step)• Diﬀusive volume flux in the 2-D continuity equation (both barotropic and baroclinic time step)• Tracer advection flux by Second Order Moment and MPDATA schemes• Tracer diﬀusion flux by horizontal diﬀusion• Tracer flux by isoneutral transport scheme (isoneutral diﬀusion plus Gent-McWilliams parameterization)• Sea ice advection plus diﬀusion flux (area fraction, volume, thermal energy) with some special treatment. See
Chapter 17 for details.

18.4.2 Surface fluxes
In order to impose conservation on some properties in a set of nested models, surface fluxes must be also taken into
account.
In MRI.COM, you may impose conservation condition for the following cases by adjusting the net surface flux.
• Volume of sea water in an ocean-only mode (WADJ option).
– 232 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 18

NestingΩωΓFigure18.6 Schematic for considering conservation in a system of coarse and fine grid models. The coarse grid(parent) model region is denoted by Ω and the fine grid (child) model region is denoted by ω. The boundary between
models is denoted by Γ. These notations follow Figure 1 of Debreu et al. (2012).

(a)

(b)

ic

ip

(c)

+

+

diﬀerence divided by+total cross sectional area

Figure18.7 Schematic explaining flux adjustment for a child model. In this zonal-vertical section, tracer points
coincide at i p and i c . (a) Meridional fluxes are computed in the child model using the values interpolated from theparent model. Flux is integrated over the corresponding parent grid and this is compared with that of the parent model(b). (c) Diﬀerence between the integrated fluxes is divided by the integration area (unit parent grid area) to give flux
correction for the child model (pink).

• Salt content in the ocean - sea ice system in any mode (l_sss_rst_cnsv = .true. in namelist nml_sss_restore)• Heat content and volume of sea water in a coupled mode (SCUPCGCM and NPUTFLUX options for parent and NGETFLUX
option for child)
For example, to conserve sea water volume in the system (WADJ option), a globally constant adjustment factor f w adj isused so that∫∫p
( f w + f w adj )dS +( f wc + f w adj )dS = 0,
(18.20)ωΩ−ωpwhere f w and f wc represent total surface water fluxes of parent and child models, respectively. The adjustment factor isgiven byf w adj = −

(

∫

∫
p

f w dS +Ω−ωω)
f wc dS /(SΩ−ω + Sω ),(18.21)∫∫
where SΩ−ω = Ω−ω dS and Sω = ω dS.
To avoid a drift of total salt content of the system (l_sss_rst_cnsv = .true. in namelist nml_sss_restore) in thecase of surface salinity restoring, a globally constant adjustment factor f s adj is used so that∫∫p
( f s + f s adj )dS +( f sc + f s adj )dS = 0,
(18.22)Ω−ωω– 233 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201718.4 Conservation
parent

child

adv

jp jc
diﬀ

×：tracer●：velocity

Figure18.8 Computation of parent model’s zonal diﬀusive flux across the boundary, coinciding at j p and jc in termsof the tracer point. In this case, the sum of diﬀusive fluxes of the child model (green) is used as the diﬀusive flux forthe parent model at j p . This is contrasted with the advective flux (blue) where the parent model flux is used to adjust
the child model flux.

pwhere f s and f sc represent surface salinity flux due to restoring in parent and child models, respectively. The adjustment
factor is given by
∫
∫
(
)
p
f s adj = −f s dS +f sc dS /(SΩ−ω + Sω ).
(18.23)ωΩ−ω

When a set of nested grid models is coupled with an atmospheric model (SCUPCGCM and NPUTFLUX options for parent
and NGETFLUX option for child), the sum of sea surface heat flux of the ocean model must equal that of the atmosphericmodel. To achieve this, a globally constant adjustment factor f q adj is used so that
∫

∫
Ω−ω

p
( f q + f q adj )dS +∫ω( f qc + f q adj )dS =f qa dS,(18.24)Ωpwhere f q , f qc , and f qa represent net surface heat flux in parent oceanic, child oceanic, and atmospheric models, respectively.
The adjustment factor is given by
∫
[∫
(∫
)]
p
f q adj =f qa dS −f q dS +f qc dS /(SΩ−ω + Sω ).
(18.25)ΩωΩ−ω

f q adj is used to oﬀset the net longwave flux in the ocean models.
These adjustments are applicable only to on-line two-way nesting by using pre-communicator of Scup.

18.4.3 Nudging
To conserve the total amount of tracer of the parent model when a tracer of the parent model is nudged toward that of theϕchild model, an additional flux adjustment is needed. A globally constant adjustment factor Sadj is added to the right handside of (18.9) so that∫1p
ϕ(18.26)[− (ϕ p − ϕC ) + Sadj ]dV = 0.
γωThe adjustment factor is given byϕ

Sadj =

∫

1
ω γ

p

(ϕ p − ϕC )dV /Vω,

– 234 –

(18.27)

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 18

Nesting∫
where Vω = ω dV . By doing so, nudging only changes the distribution of the tracer in the parent model withoutaﬀecting its total amount. MRI.COM oﬀers this adjustment as an option (flg_feedback_cnsv_ndg = .true. in namelist
nml_nest_par).

18.5 Sponge region
Appearance of noises and discontinuity around the parent-child boundary is almost unavoidable. It is a common practice
to set a sponge region around the boundary to reduce them. Several methods are available for MRI.COM.

18.5.1 Giving 2-D distribution of diﬀusion and viscosity coeﬃcient
A simple choice would be to give 2-D distribution of time-invariant coeﬃcient for diﬀusion and viscosity operators thathave been selected for each model. For example, when a common operator is selected for both parent and child models,the coeﬃcients of the child model should be the same as those of the parent model around the side boundary to obtain
similar results there. The coeﬃcients in the boundary region should be connected smoothly with the small coeﬃcients in
the interior.
Among the mixing operators available in MRI.COM, the additional harmonic viscosity is most eﬀective for getting
a smooth solution around the parent/child boundary. However, this scheme should be used with caution because theharmonic viscosity may even damp resolved processes of the parent model. These additional diﬀusivity and viscosity areavailable, whether nesting is one-way or two-way.

18.5.2 Scheme specialized for two-way nesting
MRI.COM has some schemes specialized for two-way nesting in order to make the model stable.
a. Parent model

A parent model can use an additional Laplacian operator for both velocity and tracer. This operator may be applied to thevelocity points on the boundary (red line in Figure 18.9) and to the tracer point (blue line) just outside the feedback region
where parent-model values are upscaled (replaced) by child-model values (light blue region). Laplacian is calculated in
the direction perpendicular to the boundary.

×：tracer●：velocity

Figure18.9 Same as Figure 18.1b except for indicating of sponge regions for a parent model. Laplacian diﬀusionand viscosity operators are applied in the blue and red bands, respectively.

– 235 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201718.6 Usageb. Child model

For a child model, Laplacian operator is applied to the diﬀerence between the child value and the parent value mapped on
the child grid in the sponge region following Section 2.5 of Debreu et al. (2012),∂ϕc= ∇ · κ∇(ϕc − ϕcP ),∂t(18.28)where ϕc is the child model value, ϕcP is the parent value mapped on the child grid, and κ is viscosity/diﬀusion coeﬃcient.
The operator is prepared for velocity and tracers. The width of the sponge region is two parent velocity cells from theboundary constructed by connecting tracer cells (region shaded with purple in Figure 18.10).

×：tracer●：velocity

Figure18.10 Same as Figure 18.1b except for indicating of sponge region for child model. Laplacian diﬀusion/viscosity operator is applied to the diﬀerence between child and parent values in the region shaded with purple.

18.6 Usage
18.6.1 Oﬄine one-way nesting (OFFNESTPAR, OFFNESTSUB option)a. Compilation

The model option for the parent model is OFFNESTPAR (OFFline-NESTing PARent), and that for the child model is
OFFNESTSUB. These model options and the model name should be specified in configure.in.
b. Runtime specification for the parent modelAt run time, namelists nml_parentmodel_grid (Table 18.2) and nml_parentmodel_package (Table 18.3) must be
specified. They are read by parent/parent_common.F90 on execution. See docs/README.Namelist for details.
Following files are output from the parent model to be read by the child model.
• Three rows or columns of data for each side boundary of specified packages. For package (package_name),’uv’, ’barotropic’, ’active_tracer’, ’passive_tracer’, ’vmixcoef’, ’nohkim’, ’my25’, ’gls’, ’ice_thermodynamics’,’ice_dynamics’, and ’tide’ are available.
• Land-sea index of the parent model (file_parent_sea_index_out).
• Latitude and longitude of the parent model (file_parent_grid_out).
The daily output of side boundary data will work in running the child model.

– 236 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 18

Nesting
Table18.2 Namelist nml_parentmodel_grid for the oﬀ-line parent model

variable namesub_region_ifirst

dx_ratio_w

description
the western end point in the parent model to be
used to interpolate data (T-point)
the eastern end point in the parent model to be
used to interpolate data (T-point)
the southern end point in the parent model to be
used to interpolate data (T-point)
the northern end point in the parent model to be
used to interpolate data (T-point)
refinement ratio at west boundary

dx_ratio_erefinement ratio at east boundary

dy_ratio_srefinement ratio at south boundary

dy_ratio_nrefinement ratio at north boundary

file_parent_sea_index_out
file_parent_grid_out

land-sea index of the parent model
latitude-longitude of the parent model

sub_region_ilast
sub_region_jfirst
sub_region_jlast

usage

odd number, zero for no nestingalong this side
odd number, zero for no nestingalong this side
odd number, zero for no nestingalong this side
odd number, zero for no nestingalong this side

Table18.3 Namelist nml_parentmodel_package for the oﬀ-line parent model

variable namepackage_namefile_root_namebndr_dt_sec
nrec_firstdescription
root name of the output boundary data files
time interval [sec] of output
the first record number of the output files

l_write_initial_state

usage

must be multiple of model time step
= 1 : create new file
> 1 : append to old file.true. : output initial state.false. : do not output initial statec. Runtime specification for the child modelAt run time, namelists nml_submodel_grid (Table 18.4), nml_submodel_bnd_cnd (Table 18.5), and nml_submodel_package(Table 18.6) must be specified. They are read by sub/sub_common.F90 on execution. See docs/README.Namelist for
details.
To run the child model, prepare the following two data files:• file_sub_wgt_t: The contribution ratio between parent and child models around the boundary (T-point),• file_sub_wgt_u: The contribution ratio between parent and child models around the boundary (U-point).
In addition, produce the following files by running the parent model:• file_root_name in nml_submodel_package: Three rows or columns of the necessary elements at each sideboundary,• file_parent_sea_index: Land-sea index of the parent model,• file_parent_grid: Latitude and longitude of the parent model.
Note the following when preparing barotropic boundary data. Since time filtering is used to feedback the result of thebarotropic equations to the baroclinic mode, the barotropic equations are integrated over the end of the baroclinic time.
Thus, the child model needs barotropic data longer than the original output from the parent model. To fulfill this need,prepare the barotropic data file to include one additional time record. For example, in a historical (yearly) run, the second
record of the following year of the parent model should be appended to the present year data. If the additional record is
not prepared, MRI.COM integrates the barotropic equation by repeating the final record just as needed.

– 237 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201718.6 Usage
Table18.4 Namelist nml_submodel_grid for the oﬀ-line child model

variable nameparent_grid_im

dx_ratio_w

description
Total number of grid points in X directionof the parent model
Total number of grid points in Y directionof the parent model
the western end point in the parent model
to be used to interpolate data (T-point)
the eastern end point in the parent model tobe used to interpolate data (T-point)
the southern end point in the parent model
to be used to interpolate data (T-point)
the northern end point in the parent model
to be used to interpolate data (T-point)
refinement ratio at west boundary

dx_ratio_erefinement ratio at east boundary

dy_ratio_srefinement ratio at south boundary

dy_ratio_nrefinement ratio at north boundary

file_parent_grid
file_parent_sea_index

latitude-longitude of the parent model
land-sea index of the parent model

parent_grid_jm
sub_region_ifirst
sub_region_ilast
sub_region_jfirst
sub_region_jlast

usage

odd number, zero for no nesting along thisside
odd number, zero for no nesting along thisside
odd number, zero for no nesting along thisside
odd number, zero for no nesting along thisside

Table18.5 Namelist nml_submodel_bnd_cnd for the oﬀ-line child model

variable namefile_sub_wgt_t
file_sub_wgt_u

description
filename of weighing factor for parent/child model aroundthe boundary (T-points)
filename of weighing factor for parent/child model aroundthe boundary (U-points)

usage

Table18.6 Namelist nml_submodel_package for the oﬀ-line child model

variable namepackage_namefile_root_namebndr_dt_sec
bndr_first_date
num_bndr_record

description

usage

root name of the output boundary data files
time interval [sec] of the input data
integer array(6) indicating date and time (YMDMHS) of
the first record
the last record number of the input files

18.6.2 On-line nesting (SCUPNEST option)a. Compilation

For on-line mode, SCUPNEST option should be specified. Also, the name of the model should be explicitly specified suchas NAME_MODEL = modelname. The model option for the parent model is PARENT, and that for the child model is SUB.
In ocean-only mode, surface fluxes are calculated in each model, preferably on the basis of common atmospheric state.
When the parent model is coupled with an atmospheric model (SCUPCGCM), the parent model receives surface fluxes
from the atmospheric model. The surface flux may be sent from parent to child model. In the parent (child) model,
NPUTFLUX (NGETFLUX) option must be selected.
These model options and the model name should be specified in configure.in.

– 238 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 18

Nestingb. Runtime specification: NAMELIST_SCUP

In on-line mode, parent and child models are run at the same time and Scup (simple coupler) by Yoshimura and Yukimoto(2008) is used to exchange data. User should tell the coupler how data are exchanged between parent and child models
via a namelist file NAMELIST_SCUP. A template of NAMELIST_SCUP is available from MRI.COM execution environment(MXE). Mapping tables for data transfer used by Scup can be created in MXE. User may comment out unnecessary dataexchange according to the model and runtime options chosen for a specific experiment. An example of NAMELIST_SCUP
is listed in the following. In this example, the model name of the parent model is GLOBAL and that of the child model is
NP01.


An example of NAMELIST_SCUP# pre-communicator## (parent -> child)&nam_scup_pre model_put=’GLOBAL’, model_get=’NP01’, type=’REAL8’ /&nam_scup_pre var_put=’ALONTC’, var_get=’ALONTC’, dst_get=’ALL’ /&nam_scup_pre var_put=’ALATTC’, var_get=’ALATTC’, dst_get=’ALL’ /&nam_scup_pre var_put=’ALONUC’, var_get=’ALONUC’, dst_get=’ALL’ /&nam_scup_pre var_put=’ALATUC’, var_get=’ALATUC’, dst_get=’ALL’ /# activate if l_sss_rst_cnsv = .true.
&nam_scup_pre var_put=’SFLXADJ’, var_get=’SFLXADJ’, dst_get=’ALL’ /## (child -> parent)&nam_scup_pre model_put=’NP01’, model_get=’GLOBAL’, type=’REAL8’ /# activate if l_sss_rst_cnsv = .true. and flg_trcflux_cnsv=.true.,&nam_scup_pre var_put=’SFLXSUB’, var_get=’SFLXSUB’, dst_get=’ALL’ /## main-communicator## (parent -> child) side boundary prognosticated variables (ALWAYS activate)&nam_scup model_put=’GLOBAL’, model_get=’NP01’ ,comp_put =’BAROCLIP_GLOBAL’ , comp_get =’BAROCLIS_NP01’ ,grid_put =’BCLI3DP_GLOBAL_T’, grid_get =’BCLI3DS_NP01_T’ ,fl_remap=’../data-np/rmp_tw_p2ct_3d.d’ /&nam_scup
var_put=’TRC_SIDE01’, var_get=’TRC_SIDE01’, intvl=-1, lag=0, flag=’SNP’/&nam_scup
var_put=’TRC_SIDE02’, var_get=’TRC_SIDE02’, intvl=-1, lag=0, flag=’SNP’/## (child -> parent) for replace (activate IF flg_feedback_replace=.true.)&nam_scup model_put=’NP01’, model_get=’GLOBAL’ ,comp_put =’BAROCLIS_NP01’, comp_get =’BAROCLIP_GLOBAL’,
grid_put =’BC3DS_NP01_TC’, grid_get =’BC3DP_GLOBAL_TA’,
fl_remap=’../data-np/rmp_c2pt_merge_3d.d’ /&nam_scup
var_put=’TRC_REPL01’, var_get=’TRC_REPL01’, intvl=-1, lag=0, flag=’SNP’/&nam_scup
var_put=’TRC_REPL02’, var_get=’TRC_REPL02’, intvl=-1, lag=0, flag=’SNP’/

c. Runtime specification for parent modelAt run time, namelist nml_nest_par (Table 18.7) must be specified. They are read by nest_scup_par.F90 on execution.
See docs/README.Namelist for details. If sponge layer is used for the parent model (flg_sponge = .true.), the twodimensional distribution of the diﬀusion and viscosity coeﬃcient in the sponge layer is read from a file whose name isgiven by namelist nml_par_diff (Table 18.8) and nml_par_visc (Table 18.9).

– 239 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201718.6 Usage
Table18.7 Namelist nml_nest_par for the on-line parent model

variable nameflg_carbon
flg_verbose
flg_sponge
flg_sponge_sub
flg_feedback_nudge
flg_feedback_replace
flg_volflux_cnsv

units
logicallogicallogicallogicallogicallogicallogicalflg_iceflux_cnsv

logicalflg_send_sshhdif_c2p
flg_trcflux_cnsv

logicallogicalflg_send_trchdif_p2clogicalflg_send_trchdif_c2plogicalflg_send_trcidif_p2clogicalflg_send_trcidif_c2plogicalchfbc_sub_day
scup_norecv
nest_depthday
1integer

num_sub_models
file_upscale

integercharacterdescription
Send surface CO2 flux or not
Output details of the processing
Sponge region is set just outside the boundary
Send data to child model’s sponge region
Use nudging of tracers as two-way feedback
Use replacing of all fields as two-way feedback
Impose conservation of volume flux at u(v)star
for child model
Impose conservation of sea ice flux at u(v)star
for child model
Send SSH diﬀusion flux from child to parent
Impose conservation of tracer flux at u(v)star
for child model
Send tracer lateral diﬀusion flux from parent tochild
Send tracer lateral diﬀusion flux from child toparentSend tracer isopycnal diﬀusion flux from parent
to parent
Send tracer isopycnal diﬀusion flux from childto parent
restoring time of received tracer in day
undefined value of received tracer given by scup
nesting is operated from surface to nest_depth
level
the number of child (sub) models
Name of the file that contains mask of the upscale region

usage

default = .false.
default = .false.
default = .false.
default = .false.
default = .false.
default = .false.
default = .false., choose
either p2c or c2p or neither
default = .false., choose
either p2c or c2p or neither
default = .false., choose
either p2c or c2p or neither
default = .false., choose
either p2c or c2p or neither

default = 1

Table18.8 Namelist nml_par_diff for the on-line parent model

variable namefile_sponge_par_diff

description
2D distribution of horizontal diﬀusivity justoutside the boundary

usage

Table18.9 Namelist nml_par_visc for the on-line parent model

variable namefile_sponge_par_visc

description
2D distribution of horizontal viscocity just outside the boundary

usage

d. Runtime specification for child modelAt run time, namelist nml_submodel_grid (Table 18.10) and nml_nest_sub (Table 18.11) must be specified. They areread by nest_scup_sub.F90 on execution. See docs/README.Namelist for details.
The two dimensional distribution of the diﬀusion and viscosity coeﬃcient in the sponge layer is optionally
read from a file whose name is given by namelist nml_sub_diff (Table 18.12), nml_sub_visc (Table 18.13) (ifflg_sponge_sub = .true.), nml_sub_isopyc (Table 18.14), and nml_sub_ssh_diff (Table 18.15).

– 240 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 18

Nesting
Table18.10 Namelist nml_submodel_grid for the on-line child model

variable nameparent_grid_im

dx_ratio_w

description
Total number of grid points in X direction of theparent model
Total number of grid points in Y direction of theparent model
the western end point in the parent model to be
used to interpolate data (T-point)
the eastern end point in the parent model to be used
to interpolate data (T-point)
the southern end point in the parent model to be
used to interpolate data (T-point)
the northern end point in the parent model to be
used to interpolate data (T-point)
refinement ratio at west boundary

dx_ratio_erefinement ratio at east boundary

dy_ratio_srefinement ratio at south boundary

dy_ratio_nrefinement ratio at north boundary

parent_grid_jm
sub_region_ifirst
sub_region_ilast
sub_region_jfirst
sub_region_jlast

l_parent_tripolarusage

odd number, zero for no nesting along
this side
odd number, zero for no nesting along
this side
odd number, zero for no nesting along
this side
odd number, zero for no nesting along
this side.true. : The parent model uses
TRIPOLAR grid

Table18.11 Namelist nml_nest_sub for the on-line child model

variable namefile_sub_wgt_t

unitscharacter

file_sub_wgt_ucharacter

file_sub_wgt_x

character

file_sub_wgt_ycharacter

flg_carbon
flg_verbose
flg_sponge_sub
flg_feedback_nudge
flg_feedback_replace
flg_volflux_cnsv

logicallogicallogicallogicallogicallogicalflg_iceflux_cnsv

logicalflg_send_sshhdif_c2p
flg_trcflux_cnsv

logicallogicaldescription
filename of weighing factor for parent/childmodel around the boundary (T-points)
filename of weighing factor for parent/childmodel around the boundary (U-points)
filename of weighing factor for parent/childmodel around the boundary (Ustar-points)
filename of weighing factor for parent/childmodel around the boundary (Vstar-points)
Send surface CO2 flux or not
Output details of the processing
Sponge region is set just inside the boundary
Use nudging of tracers as two-way feedback
Use replacing of all fields as two-way feedback
Impose conservation of volume flux at u(v)star
for the child model
Impose conservation of sea ice flux at u(v)star
for the child model
Send SSH diﬀusion flux from child to parent
Impose conservation of tracer flux at u(v)star
for the child model

usage

default = .false.
default = .false.
default = .false.
default = .false.
default = .false.
default = .false.
Continued on next page– 241 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201718.7 Program structure

variable nameflg_send_trchdif_p2c
flg_send_trchdif_c2p
flg_send_trcidif_p2c
flg_send_trcidif_c2p
nest_depth

Table 18.11 – continued from previous page
units
description
logicalSend tracer lateral diﬀusion flux from parent tochildlogicalSend tracer lateral diﬀusion flux from child toparentlogicalSend tracer isopycnal diﬀusion flux from parent
to parent
logicalSend tracer isopycnal diﬀusion flux from childto parent
integer
Nesting is operated from surface to nest_depth
level

usage
default = .false., choose
either p2c or c2p or neither
default = .false., choose
either p2c or c2p or neither
default = .false., choose
either p2c or c2p or neither
default = .false., choose
either p2c or c2p or neither

Table18.12 Namelist nml_sub_diff for the on-line child model (tracer_nest.F90)

variable namefile_sponge_sub_diff

description
3D distribution of horizontal diﬀusivity justinside the boundary

usage
valid if flg_sponge_sub = .true.

Table18.13 Namelist nml_sub_visc for the on-line child model (clinic_nest.F90)

variable namefile_sponge_sub_visc

description
3D distribution of horizontal diﬀusivity justinside the boundary

usage
valid if flg_sponge_sub = .true.

Table18.14 Namelist nml_sub_isopyc for the on-line child model (ipcoef.F90)

variable namel_isopyc_mask_file
file_isopyc_mask_sub

description.true./.false 2D distribution of isopycnal diﬀusion mask is read from file
2D mask of isopycnal and GM diﬀusionusage

diﬀusion coeﬃcients are determined as in the normal case

Table18.15 Namelist nml_sub_ssh_diff for the on-line child model (surface_integ.F90)

variable namel_ssh_diff_file
file_ssh_diff_sub

description.true./.false 2D-distribution of SSH diﬀusion isread from file
2D distribution of horizontal diﬀusivity of SSH

18.7 Program structure
The program structure of an oﬀ-line parent model is as follows.
ogcm__prepare|+- parent_ctl__ini|+- parent_ctl__prepare

ogcm__run|+- parent_ctl__write|+- parent_ctl__fin_run(initialize the parent model)(output boundary data of the initial state)(output boundary data at specified interval)(close files)– 242 –usage

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 18

Nesting

The program structure of an oﬀ-line child model is as follows.
ogcm__ini|+- sub_ctl__ini|+- sub_ctl__prepare(initialize the child model)(open boundary files, etc.)

ogcm__run|+- part_1|||+- surface_ctl__predict|||+- surface_predict
|
|
|
|
|
+- surfce_integ__predict
|
|
|
|
|
+- DO LOOP_N
|
|
|
|
|
|
|
+- surface_sub__input (replace barotropic boundary data)
|
|
|
|
|
+- surface_sub__input (replace time-filtered data)|||+- surface_bcli_predict_h (l_global_local_cnsv_ssh = .true.)|||+- surface_sub__input (replace time-filtered data again)|+- part_2|+- sub_ctl__input_side(read/receive baroclinic boundary data|
after "tracer_ctl__implicit_z" and before "strat_adjust_ctl__main" )|+- sub_ctl__input_side_vd (read/receive viscosity/diffusivity data
after mixed layer model)– 243 –