TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Part VI

Miscellaneous

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 19

Basics of the finite diﬀerence method
This chapter describes the basics of finite diﬀerence methods for solving diﬀerential equations. The general principles ofthe finite diﬀerencing methods are introduced using the diﬀusion equation as an example in Section 19.1. Sections 19.2 and19.3 describe applying finite diﬀerence methods of time and space derivatives in diﬀerential equations. Considerations infinite-diﬀerence methods for advection-diﬀusion equations are discussed in Section 19.4. An implicit method for solving
the diﬀusion equation is described in Section 19.5. Section 19.6 summarizes discretization schemes used by MRI.COM.
Durran (1999) treats the basics of finite diﬀerence methods for solving diﬀerential equations of advection and diﬀusionin geophysical fluid dynamics.

19.1 Diﬀusion equation
As an example, consider an initial-boundary value problem expressed by a one-dimensional diﬀusion equation (heatconductive equation),∂ 2T∂T= κ 2.
(19.1)∂t∂x
Given T (x, 0) = f (x) as the initial distribution and T (0, t) = T (L, t) = 0 as the boundary condition, the analytical solutionis∞∑2
T (x, t) =f m e−κkm t sin(k m x),(19.2)
m=0

where
fm =

2
L

∫ L
f (x) sin(k m x)dx,

km =

0

mπ
.
L

(19.3)

Next, consider the finite diﬀerence method to get the solution numerically. In the finite diﬀerence method, grids are set
with a finite increment in space and time, and each term in the equation is evaluated at each grid using Tjn = T (x j , t n ). For
example,nn
Tj+1− 2Tjn + Tj−1
Tjn+1 − Tjn=κ,(19.4)∆t∆x 2
where ∆t = t n+1 − t n and ∆x = x j+1 − x j .
Distribution at the new time level T n+1 can then be calculated if T n is known. This finite diﬀerence equation is identicalto the original diﬀerential equation (19.1) in the limit ∆t → 0, ∆x → 0 (consistency).
If the initial distribution is assumed to be f (x) = T0 sin k1 x, the solution of the finite diﬀerence equation (19.4) fort = t n is
Tjn = λ nT0 sin k1 x j ,(19.5)
where
λ =1−

2κ∆t
∆x 2

(1 − cos k1 ∆x).

(19.6)In order to suppress oscillation and divergence of the solution (stability), 0 < λ < 1 is necessary and ∆x and ∆tmust be set to satisfy this condition. This solution is identical to the analytical solution in the limit of ∆t → 0, ∆x → 0(convergence).
To summarize, the finite diﬀerence method that satisfies consistency, stability, and convergence is the necessary condition
for an accurate solution.

– 247 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201719.2 Finite diﬀerence expressions for time derivatives

19.2 Finite diﬀerence expressions for time derivatives
The following four finite diﬀerence expressions are employed for the time derivatives in MRI.COM:
T n+1 − T n= F (T n )∆t
T n+1 − T nbackward := F (T n+1 )∆t
T n+1 − T nT ∗n+1 − T n= F (T n ),= F (T ∗n+1 )Matsuno :∆t∆t
T n+1 − T n−1
leap-frog := F (T n ).
2∆tforward :(19.7)(19.8)(19.9)(19.10)The scheme used in the previous section is the forward scheme. The forward, backward, and Matsuno schemes use thevalues at two time levels and are accurate to O(∆t), while the leap-frog scheme uses three time levels and is accurate toO(∆t 2 ). Basically, the leap-frog scheme is employed in MRI.COM because of its higher order accuracy.
However, the leap-frog scheme cannot be applied to the diﬀusion equation. A solution by the finite diﬀerence method
using the leap-frog scheme is given by
Tjn = (Ta λ na + Tb λ bn ) sin k1 x j ,(19.11)√−α + α 2 + 4
λa = −
,
2

where

√
−α − α 2 + 4
4κ∆t(1 − cos k1 ∆x)).
λb = −(α ≡2∆x 2(19.12)

Because λ b < −1 for arbitrary values of α, the divergent mode with oscillation is always included (computational mode).
In order to avoid this computational mode, the forward scheme is employed for diﬀusion and viscosity terms in MRI.COM.
When the diﬀusion and viscosity coeﬃcients are very large as in the surface mixed layer, the time step has to be unusuallysmall for the stability of the forward scheme according to (19.6). In such a case, the backward scheme is used for vertical
diﬀusion and viscosity (implicit method; see Section 19.5). Though the time integration at each point can proceed without
referring to the result of other points by the forward, leap-frog, and Matsuno schemes, it must be done by solving combined
linear equations in the backward scheme (see Section 19.5).
The Matsuno scheme is useful for suppressing the computational mode in the leap-frog scheme. By defaults, theMatsuno scheme is used once per twelve steps of the leap-frog scheme in MRI.COM. This interval can be changed at run
time using a namelist parameter (nstep_matsuno_interval) of namelist nml_time_step (Table 21.6). It should be
noted that the Matsuno scheme needs twice as many numerical operations as the forward and leap-frog schemes.

19.3 Finite diﬀerence expression for space derivatives
Let us consider a one-dimensional advection equation,∂T∂T= −u ,∂t∂x(19.13)T (x, t) = T (x − ut, 0).

(19.14)

where u is a constant velocity. The solution is

Using the leap-frog scheme for time diﬀerencing, the finite diﬀerence equation can be written as follows:
Tjn+1 − Tjn−1
2∆t

= −u

n
n
Tj+
1 −T
j− 1
2

2

∆x

,

(19.15)where T n 1 and T n 1 are the values at the left and right faces of the grid cell for x j , i.e., values at x j− 1 and x j+ 1 . The point
j− 2

j+ 2

2

2

x j− 1 is defined as the central point between x j and x j−1 . Because the transport of T at the boundary that enters a grid2cell is identical to that leaving the adjacent grid cell, the total T in the whole system is conserved in this finite diﬀerence
equation.
– 248 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 19

Basics of the finite diﬀerence method

There are several methods to decide T n 1 using a value at a single or multiple grid points. The following are two simple
and frequently used formulations,

j− 2

n
n
n
n
Tj−
upstream finite diﬀerence : Tj−
1 = T j (u < 0),
1 = T j−1 (u > 0),2(19.16)2ncentral finite diﬀerence : Tj−
1 =

n
+ Tjn
Tj−1

2

2

.

(19.17)

The former is accurate to O(∆x), and the latter is accurate to O(∆x 2 ).
In central finite diﬀerencing, the expression for (19.15) is
Tjn+1 − Tjn−1

= −u

2∆t

n
n
− Tj−1
Tj+1

2∆x

.

(19.18)Assuming the solution to be T (x, t) = τ(t)e−ik x ,τ n+1 = τ n−1 + 2iατ n, where α ≡u∆tsin k∆x.
∆x(19.19)It is stable (neutral) if |α| ≤ 1. To be stable for any wave number,
u∆t
≤1
∆x

(19.20)must be satisfied (CFL condition). However, if τ n = τ 0 e−in∆θ ,∆θ = − sin−1 [µ sin k∆x](where µ ≡

u∆t
).
∆x

(19.21)

Expanding the r.h.s. by a Taylor expansion we obtain1(µ sin k∆x) 36µ(k∆x) 3 µ3 (k∆x) 3≃ − µk∆x +
−
6
6
{
}
2
(k∆x)
= − µk∆x 1 −(1 − µ2 ) .
6∆θ ≃ − µ sin k∆x −(19.22)This means that the phase of the solution from this finite diﬀerence scheme is delayed relative to that of analytical solution,
depending on its wavenumber (numerical dispersion). Therefore, a distribution with maxima and minima that do notexist in the initial distribution arises. However, this method is popularly used since the kinetic energy is conserved by
employing the central diﬀerence in the advection term in the equation of motion. Moreover, the "Arakawa method," which
can nearly conserve the enstrophy (squared vorticity) for horizontally non-divergent flows, is adopted in MRI.COM byusing the central diﬀerence. This topic is treated in Chapter 7.
Using the upstream finite diﬀerence, the finite diﬀerence equation (19.15) is
Tjn+1 − Tjn−1−un
Tjn − Tj−1.

(19.23)∂T u∆x ∂ 2T+
+ O(∆x 2 ).
∂x
2 ∂ x2(19.24)2∆tExpanding the r.h.s. by a Taylor expansion we obtain

= −u

∆x

The second term has the diﬀusion (heat conductive) form (which disappears in the central finite diﬀerencing). Actually, theinitial distribution diﬀuses when the advection equation is solved by the upstream finite diﬀerence (numerical diﬀusion).
The third order schemes (QUICK, QUICKEST, and UTOPIA) can be used in MRI.COM to suppress the numericaldispersion and diﬀusion somewhat in the advection calculation for tracers, but not completely. The grid boundary value
is set in QUICK as
n
Tj−
1 =

n
n
−Tj−2
+ 6Tj−1+ 3Tjnnn
3Tj−1+ 6Tjn − Tj+1

(u < 0).
(19.25)8
8
2The QUICKEST method uses the time averaged value at the grid boundary as the tracer value to be transported, considering
the change of the value there by advection during one time step. UTOPIA is a multi-dimensional extension of QUICKEST.
The details of these schemes are described in Chapter 8.
2n(u > 0), Tj−
1 =– 249 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201719.4 Finite diﬀerencing of advection-diﬀusion equation

19.4 Finite diﬀerencing of advection-diﬀusion equation
According to the above restriction, when the advection-diﬀusion equation (Eqs. (8.1) and (9.1)) is expressed in finite
diﬀerence form using the leap-frog scheme, it is necessary to use present (previous) time level for advection (diﬀusion)
term. The following finite diﬀerence equation is then employed:
T n+1 − T n−1= −A(T n ) + D(T n−1 ),2∆t(19.26)

where A and D are advection and diﬀusion operators, respectively.
When the vertical diﬀusion term is very large, D(T n+1 ) is used instead of D(T n−1 ). This formula is an implicit scheme
and is described in the next section.

19.5 Implicit method for vertical diﬀusion equation
Turbulent mixing is parameterized by using high vertical diﬀusivity and viscosity determined by boundary layer models,
which was treated in Chapter 15. The time step must be set very small to keep the calculation stable when viscosity anddiﬀusivity are very high, since the time tendency becomes very large due to the rapid mixing. To avoid this problem, theimplicit method uses the advanced (mixed) state for evaluating viscosity and diﬀusivity, unlike the normal explicit method
where previous or present values are used.
Expressing the present time step as n and the time step before and after as n ± 1, the finite-diﬀerence method is applied
to the advection-diﬀusion equation using the leap-frog scheme. The diﬀusion term is written separately using (n − 1) stepfor the horizontal direction and (n + 1) step for the vertical direction,(T n+1 − T n−1 )= −A(T n ) + D H (T n−1 ) + DV (T n+1 ).
2∆t(19.27)

Putting all the terms involving T n+1 on the l.h.s.,
T n+1 − 2∆tDV (T n+1 ) = T n−1 + 2∆t(−A(T n ) + D H (T n−1 ))(19.28)

is obtained, which is an algebraic equation for T n+1 . The finite diﬀerence form is rewritten specifically as
Tkn+1 − 2∆t

)
1 (
n+1
n+1
κ k− 1 (Tk−1− Tkn+1 )/∆zk− 1 − κ k+ 1 (Tkn+1 − Tk+1)/∆zk+ 1
2
2
2
2
∆zk
n−1
n
= Tk + 2∆t(−A(Tk ) + D H (Tkn−1 )).

By puttinga=2∆t κ k− 1

2

∆zk ∆zk− 1,

b = 1 + a + c,c=22∆t κ k+ 1

2

∆zk ∆zk+ 1,(19.29)(19.30)2

we get,n+1
n+1= Tkn−1 + 2∆t(−A(Tkn ) + D H (Tkn−1 )).
+ bTkn+1 − cTk+1−aTk−1(19.31)

Setting Fk ≡ −A(Tkn ) + D H (Tkn−1 ), this is expressed in the matrix form as
b
*. −a
..
..
..
..
..
.
,

−c
b
−a

−c
b
..
.

−c
..
.

−a

..

.

b
−a

−c
b
−a

T n+1
T n−1 + 2∆tF1+/+/ *. T1n+1 +/ *. T1n−1 + 2∆tF
2/
2// .. 2n+1 // ..
// .. T3n−1 + 2∆tF3 ///// .. T3// .
..
// .. .. // = ..
. / .
.
//
// .. n+1
/ . n−1// .. TKM−2 // .. TKM−2 + 2∆tFKM−2 //−c / .. T n+1 // .. T n−1 + 2∆tFKM−1 //
KM−1

n+1
b - , TKM-

The l.h.s. has the form of the tri-diagonal matrix.

– 250 –

KM−1

n−1
, TKM + 2∆tFKM -(19.32)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 19

Basics of the finite diﬀerence method

19.5.1 A solution of tri-diagonal matrix
In general, simultaneous linear equations for n variables with tri-diagonal matrix coeﬃcients
B1
*. A
.. 2
..
..
,

C1
B2
..
.

C2
..
.
An−1

..

.

Bn−1
An

X1
D1
+/ *. X +/ *. D +/// .. .2 // .. . 2 //// .. .. // = .. .. ///.
/ .
/Cn−1 / . X n−1 / . Dn−1 /Bn - , X n - , Dn -(19.33)are solved using the Thomas method, which is modified from LU decomposition,
P1 = C1 /B1
Q1 = D1 /B1Ck
Pk =Bk − Ak Pk−1
Dk − Ak Q k−1
Qk =Bk − Ak Pk−1
Xn = Qn
Xk = Q k − Pk Xk+1

(19.34)(19.35)(2 ≤ k ≤ n − 1)

(19.36)(2 ≤ k ≤ n)

(19.37)(1 ≤ k ≤ n − 1).

(19.38)(19.39)19.6 Summary of schemes used by MRI.COM
MRI.COM employs centered diﬀerence schemes for spatial discretization except for the tracer advection terms.
Figure 19.1 shows the schematic of the relation between components of the governing equations in terms of time
integration. The situation depicted here is the advancement in time from n to n + 1, represented by the thick red arrows.
The main components employ the leap-frog scheme, inserting Matsuno (Euler-backward) scheme with a specified interval.
Sea ice and turbulence closure employs forward scheme.
Colored arrows represent data transfer among the components in advancing in time using the leap-frog time (n → n + 1).
(a) State at the starting time (n − 1) is used for computing surface fluxes. (b) Surface current at the starting time and SST
and SSS of the present time is given to the sea ice component. (c) Flux from the sea ice component is averaged withthat of the previous time step for use in the ocean component (d) for conservation. (e) Horizontal velocity field of thepresent time is used to compute vertical velocity and tracer transport. (f) Diagnosed vertical velocity is used immediately
in the advection schemes of momentum and tracers. (g) Temperature and salinity at the present time is used to computethe pressure gradient term. (h) Vertically integrated forcing term computed by the baroclinic component is used by thebarotropic component. (i) The barotropic transport is returned to the baroclinic component and used immediately as thevertically averaged velocity for computing the total velocity. For turbulent closure scheme, states at both present (j) andnew (k) time level are used to estimate shear production and buoyancy loss of turbulent kinetic energy. (l) Turbulent
kinetic energy at the present time is used to compute vertical mixing coeﬃcient. Blue arrows represent the transfer of
extant data obtained by the previous time steps. Green arrows represent the immediate use of the new data computed inthis time step. These relations constrain the order of operations within a time step. The source side of the green arrows
must be calculated earlier than the recipient side.

– 251 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201719.6 Summary of schemes used by MRI.COM

Result of this time step is used immediately
Advance of variablesat the time step nResults of previous time steps are used

n-1

n

n+1

n+2

Sea ice(c)(c)

Surface flux
(τ,Ft,Fs,Fw)

(d)

(a)

(a)

(b)

(b)

Contintuity(w)(f)(f) (e) (e)Baroclinic
momentum(u,v)
Barotropic
momentum
(U,V)

(i)

(h)

(e)

(e)

(g)

Surfaceheight(η)(j) (k)Tracers(θ,s)(j)

Turbulence
closure
(q)

(k)

(l)

Figure19.1 Schematic of the relation between components of the governing equations in terms of time integration.
The situation depicted here is the advancement in time from n to n + 1, which is represented by the thick red arrows.
Blue arrows represent the transfer of extant data obtained by the previous time steps. Green arrows represent theimmediate use of data computed at this time step. These relations constrain the order of operations within a time step.
See text for the detailed explanation.

– 252 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 20

Generalized orthogonal curvilinear coordinate
gridsThis chapter introduces generalized orthogonal curvilinear coordinates and presents related calculus.

20.1 Outline
A finite volume ocean model on geographic coordinates does not have any problem concerning the South Pole because
it does not calculate around the South Pole. However, serious problems arise around the North Pole where the meridian
concentrates to one point in the ocean. First, it is necessary to calculate the temporal evolution of the physical quantity ina special way only there, because the relations between U-cells that surround the North Pole and the northernmost T-cellare topologically peculiar. Next, even if a cell doesn’t touch the North Pole, its zonal lattice interval is extremely small
near the North Pole. Therefore, a short time step for integration is required owing to the limitation of the CFL condition.
This limitation is reflected directly in the increased calculation time required. Moreover, when the zonal grid intervals in
low latitudes and the Arctic region are extremely diﬀerent, the arguments about accuracies of numerical schemes and theparameters for diﬀusion and viscosity operators generally cannot be applied uniformly to a model domain.
The following can be considered to avoid such problems concerning the North Pole. 1) Creating a huge island including
the North Pole. The finite-diﬀerence calculation in the island is abandoned, and the lateral boundary values are restored
to the climatology. 2) Shifting the singular points of the model to a continent or a huge island by changing the model’s
horizontal grid system. The MRI.COM scheme adopts the latter approach, which is outlined in this section.
Because the MRI.COM code is based on generalized orthogonal coordinates, the geographical latitude (ϕ) and longitude(λ) are not of a great concern for the calculus in the model. However, it is necessary to know the land and sea distribution,
sea depth, scale factor, and the Coriolis parameter given as a function of λ and ϕ at every grid point of the model prior
to the calculation. We describe the method of generating the orthogonal coordinate grid system in Section 20.2. Using a
conformal transformation in the general sense, the functions that describe the relation between model coordinates (µ, ψ)and geographic coordinates (λ, ϕ), λ(µ, ψ), ϕ(µ, ψ), µ(λ, ϕ), and ψ(λ, ϕ), are obtained.
Because an atmospheric boundary condition is given in many cases at grid points in geographic coordinates, it isnecessary to prepare tables for converting the surface atmospheric temperature, the wind stress, and so on. To convert avector quantity, we must remember that the direction of the µ contour diﬀers from that of the λ contour (meridian). The
diﬀerence is described in Section 20.3. We can use the functions, (λ, ϕ) ⇐⇒ (µ, ψ), to convert a scalar quantity as well
as sea depth and the Coriolis parameter. The total flux that the ocean receives from the atmosphere should be equal to thetotal flux that the atmosphere gives to the ocean. The method for conserving the total flux is explored in Section 20.4. The
vector operation in generalized orthogonal coordinates is concisely described in Section 20.5.

20.2 Generation of orthogonal coordinate system using conformal mapping
We designate the plane that touches the sphere at the North Pole as SN . A polar stereographic projection is a conformaltransformation in the general sense, so that an orthogonal coordinate system on the sphere is mapped onto an orthogonal
coordinate system on SN and the orthogonality is preserved on the reverse transformation (Figure 20.1). Moreover, if SN isassumed to be a complex plane, various conformal transformations can be defined on it. Therefore, applying (i) the polarstereographic projection, (ii) a conformal transformation on the complex plane SN , and (iii) the reverse polar stereographic
projection to a geographic coordinate grid point (λ, ϕ) on the sphere, an orthogonal coordinate grid point on the spherecan be obtained (Bentzen et al., 1999).
The functions µ(λ, ϕ) and ψ(λ, ϕ) are obtained by the following procedure:– 253 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201720.2 Generation of orthogonal coordinate system using conformal mapping

Figure20.1 Schematic illustration of a Polar stereographic projection (a conformal transformation in the general sense
between the sphere and SN ).

1. From a point (λ, ϕ) on the sphere to a point z on SN (polar stereographic projection). Defining colatitude
ϕ ′ = π/2 − ϕ,( ′)ϕz = taneiλ,(20.1)
2
where the origin of SN corresponds to ϕ ′ = 0 (ϕ = π/2), and the positive part of the real axis corresponds to λ = 0.
2. Conformal transformation MC on SN :
ζ = MC (z).
(20.2)
3. From a point ζ on SN to a point (µ, ψ) on the sphere (reverse polar stereographic projection).
µ = arg(ζ ),(20.3)′

ψ = 2 arctan |ζ |,′ψ = π/2 − ψ .
Functions λ(µ, ψ) and ϕ(µ, ψ) are obtained by reversing the above procedure.
Defining ψ ′ = π/2 − ψ,( ′)ψ
ζ = taneiµ,2(20.4)(20.5)(20.6)

z = MC−1 (ζ ),(20.7)

λ = arg(z),(20.8)and′ϕ = 2 arctan |z|,′ϕ = π/2 − ϕ .

(20.9)(20.10)Thus, when a model coordinate grid point, (µ0 + ∆µ × (i − 1), ψ0 + ∆ψ × ( j − 1)) is given, we know the geographic positionof the point, Coriolis parameter, etc., at once.
Bentzen et al. (1999) used the linear fraction conversion as a conformal transformation on SN . That is,ζ=MC (z)=– 254 –(z − a)(b − c),(c − a)(b − z)

(20.11)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 20

Generalized orthogonal curvilinear coordinate grids

where the three complex numbers a, b, and c expressed by)
ϕ a′ iλ aa = tan
e ,
2

(

(

b = tan

ϕb′

)

2

(

e

iλ b

,

)
ϕc′ iλcc = tane ,2(20.12)

correspond to the three geographic coordinate grid points (λ a, ϕ a ), (λ b, ϕb ), and (λ c, ϕc ), which are mapped to the model
coordinate grid points, (µ, ψ) = (0, π/2), (0, −π/2), (0, 0), respectively. Therefore, the singular point (µ, ψ) = (0, π/2)
in the model calculation can be put on Greenland, by setting (λ a, ϕ a ) at 75◦ N and 40◦ W. If TRIPOLAR or JOT option is
specified instead of SPHERICAL, then two singular points: (µ, ψ) = (0, π/2) and (0, −π/2) can be put on arbitrary land
locations by suitably setting (λ a, ϕ a ) and (λ b, ϕb ), which are model parameters (north_pole_lon, north_pole_lat,
south_pole_lon, and south_pole_lat in degree) that should be specified in namelist nml_poles (Table 3.2).
When TRIPOLAR option is specified, the parameters are set to ϕ a = ϕb = 64°N, λ a = 80°E, and λ b = 100°W. The
transformed grids are used for the region north of 64°N, and geographic coordinates are used for the region south of 64°N.
This tripolar coordinate system can express the Arctic Sea with a higher resolution than the Southern Ocean. The adoptionof geographic coordinates south of 64°N enables us to do the assimilation and analysis with relative ease (Figure 20.2).
When JOT option is specified, the Joukowski conversion is used as a conformal transformation on SN . That is,ψ ′2
z = MC−1 (ζ ) = *ζ + 0 + eiµ0 .
ζ ,(20.13)This Joukowski conversion maps the area outside the circle with a radius of ψ0′ centered at the origin of the ζ-plane to thewhole domain of the z-plane, and rotates it by µ0 . The left panel of Figure 20.2 presents an example where the coordinate
system is created by setting ψ0′ to 20° and µ0 to 80°. Because there is no discontinuity of grid spacing in this coordinate
system, the singular points can be put at various positions. For instance, the singular point on the North American side
can be put on the Labrador Peninsula or in Greenland.

Figure20.2 Model coordinate grid arrangement in the Arctic sea. Left: Grid system made through the Joukowski
conversion (JOT). Right: Combination of the coordinate systems made through the linear fraction conversion andconventional geographic coordinates (TRIPOLAR).

Functions λ(µ, ψ) and ϕ(µ, ψ) are defined as subroutine mp2lp, and functions µ(λ, ϕ) and ψ(λ, ϕ) are defined as subroutine lp2mp. Module programs trnsfrm.{spherical, moebius, tripolar, jot}.F90 contain these internal subroutines.
These functions, especially mp2lp, are frequently used when the topography and the surface boundary condition are madebefore starting the main integration of model.

20.3 Rotation of vector
A vector expressed in geographic coordinates (λ, ϕ) should be rotated when observed from model coordinates (µ, ψ).
Taking advantage of the local orthogonality of coordinate axes, the angle (α) is obtained as an angle at which the meridian– 255 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201720.4 Mapping a quantity from geographic coordinates to transformed coordinates(λ = λ 0 ) of geographic coordinates intersects that of the model coordinates (µ = µ0 ) at a certain point. First, we set
z = f (ζ ),f ′ (ζ )=z = x + iy,
∂x
∂y
+i
∂u
∂u

ζ = u + iv,
=

∂y
∂x
−i .
∂v
∂v

(20.14)(20.15)

At a certain point z0 = f (ζ0 ) on geographic coordinates, the angle θ at which a curve v = v0 meets a straight line y = y0
is given by
[ / ]
[ ]
∂y ∂x
∂y
=
,
tan θ =
∂ x v0∂u ∂u v0
then (see Figure 20.3),

θ = arg( f ′ (ζ0 )).

(20.16)

Figure20.3 A meridian (red) in geographic coordinates (λ, ϕ) (left) and a meridian (blue) in model coordinates (µ, ψ) (right).

Assuming λ = arg(z) and µ = arg(ζ ), at point ζ0 the straight line (v = v0 ) meets the straight line (µ = µ0 ) at the angle−µ0 (µ0 → v0 ), and at point z0 the meridian (λ = λ 0 ) meets the curve (v = v0 ) at the angle λ 0 − θ (v0 → λ 0 ). The
meridian (λ = λ 0 ) in λ-ϕ coordinates meets the line (µ = µ0 ) in µ-ψ coordinates at angle α given as follows:
α = − µ0 + λ 0 − θ= λ 0 − µ0 − arg( f ′ (ζ0 )).

(20.17)

Subroutine rot_mp2lp defined in trnsfrm.{spherical, moebius, tripolar, jot}.F90 returns (cos α, sin α) at a
specified grid point of the model. A wind stress vector (τx, τy ) in geographic coordinates should appear in the model
ocean described in the µ-ψ coordinate system as (τx cos α − τy sin α, τx sin α + τy cos α).

20.4 Mapping a quantity from geographic coordinates to transformed coordinates
We consider a method to receive a quantity G I,J given at the geographic coordinate grids (I, J) as the quantity H M, N
at the model coordinate grids (M, N ) (Figure 20.4). The quantities are wind stress components after the vector rotation,
precipitation per unit area, sea surface atmospheric temperature, and so on. In addition, the average depth at a model grid
point can also be calculated by the following method because bottom topography (depths of sea floor) is usually given in
geographic coordinates.
– 256 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 20

Generalized orthogonal curvilinear coordinate grids

Figure20.4 Grids (I, J) and (M, N ) subdivided into finer grids (i, j) and (m, n).

Grids (I, J) and (M, N ) are suitably subdivided into finer grids (i, j) and (m, n). We call these filter grids. A quantity
Gi,′ j is assumed to be homogeneously distributed in the geographic filter grids (i, j) covered by grid (I, J),
Gi,′ j = G I,J .
′Assume the quantity at a model filter grid Hm,nis equal to that at the nearest geographic filter grid,′′Hm,n= Gi(m,n),j (m,n) .

The quantity at model grid (M, N ) is obtained as the area-weighted average:
H M, N =1∑

AH M, N m,n′
AH ′ m,n Hm,n,(20.18)where AH M, N is the area of model grid and AH ′ m,n is the area of model filter grid.
When the grid intervals of geographic filter grid (i, j) and model filter grid (m, n) are extremely small, the total quantity(flux) received on the model grids (M, N ) is equal to the total quantity (flux) given by the geographic grids (I, J). The
relation between the quantity in the geographic grids and that in the model grids is defined by weight w,∑
H M, N =w(M, N, I, J)G I,J .
(20.19)I,J

How is the quantity converted in an actual calculation in the model?
1. When the strict conservation of quantity (flux) is necessary:
Fresh water is not permitted to be generated or vanish at the surface boundary in a run using an atmosphere-ocean
coupled model, for instance. In this case, w(M, N, I, J) is prepared beforehand, and the flux is passed from theatmosphere through equation (20.19) to the ocean. The resolution of the filter grid need not be extremely fine,provided that every geographic filter grid is linked to more-than-zero model filter grids and∑∑∑∑
AG I,J =A′ G′ i, j =A′ H ′ m,n =AH M, N .
I,J

i, j

m,n

M, N

2. When conservation need not be guaranteed:
When the ocean model is driven by the surface boundary condition based on atmospheric re-analysis data, theamount of fresh water entering the sea as precipitation and river discharge is not equal to that drawn from the ocean– 257 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201720.5 Vector operation and diﬀerentiation in generalized orthogonal coordinates
through evaporation and sublimation. Therefore, the global sea surface height rises or descends during years ofintegration. It is not very important to pursue complete conservation of fresh-water under this condition. In such acase, the flux at a model grid point can be prepared beforehand using equation (20.19), to avoid the time-consuming
flux conversions in the model calculation.

20.5 Vector operation and diﬀerentiation in generalized orthogonal coordinates
To formulate the model equations, we have to know the vector operation and diﬀerentiation in generalized orthogonal
coordinates. Some basic formulae used in formulating primitive equations are presented here.
The line element vector δx at a certain point (µ, ψ, r) in an arbitrary general orthogonal coordinate system is expressedasδx = hµ δ µeµ + hψ δψeψ + hr δrer ,(20.20)where basis vectors eµ , eψ , and er are mutually orthogonal unit vectors, and hµ , hψ , and hr are scale factors.
Defining
eµ ∂
eψ ∂
er ∂
∇=
+
+
,
hµ ∂ µ hψ ∂ψ hr ∂r
the gradient of scalar A(µ, ψ, r) is∇A =eµ ∂ A eψ ∂ A er ∂ A++,
hµ ∂ µ hψ ∂ψ hr ∂r(20.21)(20.22)and the divergence of vector A = Aµ eµ + Aψ eψ + Ar er is∇·A=1h µ hψ hr

The r component of curlA is[]∂(hψ hr Aµ ) ∂(hr hµ Aψ ) ∂(hµ hψ Ar )
+
+
.
∂µ
∂ψ
∂r

[
]
∂(hψ Aψ ) ∂(hµ Aµ )1−.
h µ hψ
∂µ∂ψ(20.23)(20.24)The calculation of velocity advection includes (a · ∇)A, where a is an arbitrary vector (a = aµ eµ + aψ eψ + ar er ).
The µ component of (a · ∇)A is
(
)
(
)
Aψ
∂hµ
∂hψ
∂hµ
Ar
∂hr
a · ∇Aµ +
aµ
− aψ
+
aµ
− ar
.
(20.25)h µ hψ
∂ψ∂µhr h µ∂r∂µ
The second and third terms are so-called "metric" terms in the equation of motion in spherical coordinates.
These expressions in spherical coordinates (λ, ϕ, r) are shown next. Defining longitude λ, latitude ϕ, and radius of theearth r, scale factors are hλ = r cos ϕ, hϕ = r, and hr = 1.
Velocity vector v is
v = ueλ + veϕ + wer ,(20.26)
where eλ , eϕ , and er are the eastward, northward, and upward unit vectors, respectively, and (u, v, w) = (r cos ϕ λ̇, r ϕ̇, ṙ).
The gradient of scalar function A(λ, ϕ, r) is,∇A =∂Aeλ ∂ A e ϕ ∂ A++ er,r cos ϕ ∂λr ∂ϕ∂r

(20.27)

∇=

eϕ ∂
eλ
∂
∂
+
+ er .
r cos ϕ ∂λr ∂ϕ∂r

(20.28)where

For vector A = Aλ eλ + Aϕ eϕ + Ar er , the divergence is
∇·A=

[
]
∂ Aλ ∂(cos ϕAϕ )∂(r 2 Ar )1++.
r cos ϕ ∂λ
∂ϕr 2 ∂r

(20.29)[]∂ Aϕ ∂(cos ϕAλ )1−.
r cos ϕ ∂λ
∂ϕ(20.30)and the r component of curlA is[curlA]r =– 258 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 20

Generalized orthogonal curvilinear coordinate grids

The λ component of (a · ∇)A is[(a · ∇)A]λ = a · ∇Aλ −Aϕ aλ tan ϕ Ar aλ+.
rr(20.31)The Coriolis force in generalized orthogonal coordinates (µ, ψ, r) is given as2Ω × v = (2Ωψ w − 2Ωr v)eµ + (2Ωr u − 2Ωµ w)eψ + (2Ωµ v − 2Ωψ u)er ,(20.32)where Ω = Ωµ eµ + Ωψ eψ + Ωr er is the rotation vector of the Earth, and v = ueµ + veψ + wer is the velocity vector. We
designate f µ = 2Ωµ , f ψ = 2Ωψ , and f = f r = 2Ωr in Chapter 2. The rotation vector of the Earth is (Ωλ, Ωϕ, Ωr ) =(0, Ω cos ϕ, Ω sin ϕ) in geographic coordinates (λ, ϕ, r).

– 259 –

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 21

User’s Guide
This chapter briefly explains the procedures needed to run MRI.COM. The description in this chapter is based on MRI.COM
version 4.4 (MRICOM-4_4_beta01-20170118) and some contents presented in this chapter may not be used for the latest
version. It is recommended that users refer to README.First, README.Options, README.Namelist, README.Monitor,
and README.Restart in the docs directory when setting up a model.
The minimal information to prepare, run, and post-process is presented in this chapter in the following order:• Model setup: User defined parameter files and compilation (Section 21.1 and Table 21.1).
• Input data: Grid spacing, topography, and surface forcing etc., to be read at run time (Section 21.2 and Table
21.2).
• Restart file: Explanation of restart files (Section 21.3).
• Execution: Explanation of the runtime parameters that control time-integration (Section 21.4).
• Post process: A description of monitor files (Section 21.5).
Note that cgs units are employed to express physical values in the model.
We are developing a comprehensive package of tools "MRI.COM eXecution Environment (MXE)" that aggregatesprograms for preprocessing, execution, postprocessing, and analysis of MRI.COM experiments. This is briefly explainedin Section 21.6

21.1 Model setup
21.1.1 Model configuration file (configure.in)This section describes the procedures necessary for setting up the model and compiling its programs. First, prepare
configure.in that contains the information about model options and grid size. This is needed for compilation. Following
is an example of configure.in.
An example configure.in for Global tripolar 1° × 0.5° grid model
DEFAULT_OPTIONS="IDEALAGE ICE SIDYN CALALBSI SMAGOR VIS9P DIFAJS GLS VVDIMPSOMADVEC ISOPYCNAL HFLUX TAUBULK WFLUX RUNOFF SFLUXR LWDOWN BULKNCAR BULKITERCYCLIC ZSTAR BBL TRIPOLAR PARALLEL"#
NAME_MODEL=’GLOBAL’
IMUT=364
JMUT=368
KM=51
KSGM=51
KBBL=1
NPARTX=8
NPARTY=4
NUM_ICECAT=5
NUMTRC_P=1

– 261 –

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201721.2 Preparation of input data files for execution
Additional parameters are required for some particular model options. Those are listed in Table 21.1 (see also
README.Options).
Table21.1 Model parameters to be set in configure.in

option namealways required

variable name
IMUT, JMUT, KM
NAME_MODELKSGM
NSFMRGNBBLPARALLELKBBL
NPARTX, NPARTY

passive tracers
ICE

NUMTRC_P
NUM_ICECAT

description
zonal/meridional/vertical grid number
name of the model (default = "tmOGCM")
The number of variable layers near the sea surface. KSGM must equal KM for
ZSTAR option. See Chapter 6.
the number of side-boundary ghost cells to reduce the communication costin parallel computation (see Ishizaki and Ishikawa, 2006)
the number layers of bottom boundary layer model; must be 1
the number of zonally/meridionally partitioned region for a computation using parallel processors: the number of parallel processes should be NPARTX× NPARTY
only when any passive tracer is calculated (NUMTRC_P ≥ 1)
the number of thickness categories of sea ice

21.1.2 Compilation of the modelA standard compiling script is prepared as compile.sh in the src directory. The part depending on the system (OS,Fortran compiler, and compiler option) are found in machines directory.
To compile the programs, execute compile.sh.
The function of compile.sh is to create param.F90,
icecat/ice_param.F90, and Makefile from configure.in, param.F90.in, icecat/ice_param.F90.in, andMakefile.in respectively, by running configure and to execute the command make to create the executable file ogcm.
The environment variables for compilation are set in configure using the options prescribed in configure.in. If
configure.in is newer than param.F90, the parameter values defined in configure.in replace those in param.F90.into create new param.F90.
The program files that should be compiled are automatically selected according to the descriptions of the relationships in
Makefile, but users should be careful since it might not be perfect. The compilation should be carried out after executing./compile.sh clean when any compile option in configure.in is rewritten.

21.2 Preparation of input data files for execution
According to user’s specification of model options and job parameters, the data files listed on Table 21.2 should be preparedalong with the files that must be always prepared. See Section 21.3 for how to handle restart files.
Table21.2 Main input data files and their related program files. Here, name_model represents the specific name given
as NAME_MODEL in configure.in.

subject
runtime job parametersspecification about monitoring
variable horizontal grid spacing
variable vertical grid spacing
grid cell area and line elements
topographyreference data for tracers
restoring coeﬃcient for tracers
surface forcing
surface forcing grid
restart files

file name specified in(NAMELIST.name_model)
NAMELIST.name_model
NAMELIST.name_model.MONITOR
file_dxdy_tbox_deg (nml_horz_grid)
file_dz_cm (nml_vert_grid)
file_scale (nml_grid_scale)
file_topo (nml_topo)
trcref(_surf)_conf%file_data
rstcoef(_surf)_conf%file_data
file_data (nml_force_data)
file_data_grid (nml_force_data)
nmlrs_element, element being elements– 262 –

read fromeach package
history.F90 etc.
gridm.F90
gridm.F90
gridm.F90
topo.F90
restore_cond.F90
restore_cond.F90
force.F90
force.F90
history.F90

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

User’s Guide

21.2.1 Grid spacing and cell area
Details on how to specify grid information and how to prepare the necessary data is given in Section 3.6. The grid spacing
data file should be prepared for each of the horizontal (file_dxdy_tbox_deg) or vertical (file_dz_cm) directions when
variable grid spacing is used for that direction. The units are in degrees for the horizontal and in cm for the vertical.
When the model grid points are defined on the basis of general orthogonal coordinates, the quarter cell area and line
elements should be prepared. The units are in cgs. It is read from the file file_scale.
When spherical coordinates are used (SPHERICAL option), e.g., the grids are defined on geographical latitude andlongitude, the grid information is analytically calculated in gridm.F90, and the file file_scale is not necessary.

21.2.2 Topography
Land-sea distribution and sea-floor topography should be given by the topography data file file_topo. The topographic
data consist of the 4-byte integer array HO4(imut, jmut) that contains the sea floor depths of the velocity grid points (incm) and the 4-byte integer array EXNN(imut, jmut) that contains its corresponding vertical level. They should be written
unformatted and sequentially as follows:Format of topographic data (file_topo)
integer(4) :: ho4(imut,jmut),exnn(imut,jmut)
open(unit=nu,file=file_topo)
write(unit=nu) ho4, exnn
An example of the topography for global 1◦ × 0.5◦ model is shown in Figure 21.1. In creating a model topography,
especially for a low-resolution model, the user should be careful that the important gateways for the ocean circulation be
kept open and that the land blocking the ocean circulation be kept closed.

Figure21.1 Example of ocean model topography (global 1◦ × 0.5◦ grid model).

21.2.3 Reference data and restoring coeﬃcient for tracers
For each tracer, the integration may be started from an initial state based on the reference data and the tracer values may berestored to the reference data during the integration. To do this, tracer reference values and restoring coeﬃcients should
be prepared. See Chapter 13 for how to prepare these data.

– 263 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201721.3 Restart file

21.2.4 Surface forcing data
See Chapter 14 for how to prepare surface forcing data. The surface forcing data are read at a uniform time interval. A
leap year is set according to the calendar subroutine. Climatological data may be used repeatedly.
The following data files should be prepared according to the chosen model options. Each file is opened only once at the
beginning of run time and thus should contain all the data needed for that run.
Table21.3 Surface data to be read from namelist nml_force_data.

name

units
dyn cm−2
dyn cm−2
cm s−1cm s−1erg s−1 cm−2 =10−3 W m−2
erg s−1 cm−2

X-ward wind stressY-ward wind stress
X-ward wind speed
Y-ward wind speed
Downward shortwave radiation

U-windV-windU-windV-windShortWave

Downward longwave radiation

LongWave

Sea surface temperature
Surface air temperature
Surface air specific humidity
Scalar wind speed
Sea level pressure
precipitation
river discharge
Sea ice area fraction

SST
TempAir
SphAirScalarWind
SeaLevelPressure
Precipitation
RiverDischargeRate
IceConcentrationClimatology

◦C
◦C

1
cm s−1hPa
g s−1 cm−2
g s−1 cm−2
1

usage
if not TAUBULK
if not TAUBULK
if TAUBULK
if TAUBULKNet longwave by default(SST is required).
Downward if LWDOWN
unnecessary if LWDOWN

unnecessary if TAUBULK
also available for SLP
WFLUX
RUNOFFICECLIM

21.3 Restart file
A set of restart files provides an instantaneous state necessary to resume a model integration. This section explains how tohandle restart files. Table 21.4 lists restart variables of MRI.COM as well as their namelist that specifies their restart fileattributes.
Table21.4 Restart variables and their namelist block that specifies their restart file attributes.

variable nametemperaturesalinityx velocity
y velocity
sea surface heightx barotropic transport
y barotropic transport
x transport due to SSH diﬀusion
y transport due to SSH diﬀusion
horizontally averaged density and pressure used
for the equation of state

namelist block
nmlrs_t
nmlrs_s
nmlrs_u
nmlrs_v
nmlrs_ssh
nmlrs_uml
nmlrs_vml
nmlrs_ssh_dflx_x
nmlrs_ssh_dflx_y
nmlrs_density

model options

required if CALPP
Continued on next page

– 264 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 21

User’s Guide

Table 21.4 – continued from previous pagevariable namenamelist block
vertical diﬀusivity coeﬃcient
nmlrs_vmix_avd
vertical viscosity coeﬃcient
nmlrs_vmix_avm
passive tracers
nmlrs_ptrc
x tidal transport
nmlrs_tide_um
y tidal transport
nmlrs_tide_vm
tidal sea surface height
nmlrs_tide_ssh
vertical diﬀusion coeﬃcient of turbulence kinetic nmlrs_my_avq
energyturbulence velocity scale
nmlrs_my_q
turbulence length scale
nmlrs_my_alo
vertical diﬀusion coeﬃcient of turbulence kinetic nmlrs_nk_avq
energyturbulence kinetic energy
nmlrs_nk_eb
vertical diﬀusion coeﬃcient of turbulence kinetic nmlrs_gls_avk
energyvertical diﬀusion coeﬃcient of generic variable
nmlrs_gls_avp
turbulence kinetic energy
nmlrs_gls_eke
generic variable
nmlrs_gls_psi
turbulent length scale
nmlrs_gls_alo
moments of second order moment advection nml_somadv
schemesea ice initial state
sea ice final state

nml_seaice_rst_in
nml_seaice_rst_out

model options
VVDIMP
VVDIMPNUMTRC_P > 0
TIDE
TIDE
TIDE
MELYAMMELYAMMELYAM
NOHKIM
NOHKIM
GLS
GLS
GLS
GLS
GLS
SOMADVEC andadv_scheme%name = ”som”
in nml_tracer_data
for any tracer
ICE
ICE

21.3.1 Restart for the ocean model
For restart variables of the main part of ocean model, the attributes of their restart files are given by the namelist blocks
that start with "nmlrs_". The element names that follow are t, s, u, v, ssh, etc. as listed on Table 21.4. Requiredelements depend on model options. See docs/README.Restart for more information. Namelists must be written in
NAMELIST.name_model (default) as follows:An example namelist for temperature restart files&nmlrs_t
fname = ’result/rs_t’,/
where fname specifies the basename of the restart files. This information is shared by both input and output files. In theabove example, the name of a restart file for temperature is result/rs_t.YYYYMMDDHHMMSS. A suﬃx indicating the date
and time of data is always added to the basename.
For passive tracers, namelist nmlrs_ptrc should be repeated as many times as the number of passive tracers(=numtrc_p) with a right order. The order is determined at tracer_vars.F90. An example for a set of biogeochemical
tracers is shown in Section 11.6.
Restart files are saved in a Fortran direct-access format, which can be read by a following program.
Program to read a restart file of the ocean modelcharacter(14):: date = ’20010101000000’
real(8):: d(imut,jmut,km)
integer(4),parameter :: nu = 10 ! device number!- for 0:00z1JAN2001

open(nu, file=’result/rs_t.’//date, form=’unformatted’, && access=’direct’,recl=imut*jmut*km*8)
read(nu,rec=1) dclose(nu)– 265 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201721.3 Restart file
The file endian is that of the computer where the model runs. There are support tools to visualize restart files inthe directory, tools/MK_GRADS. Please follow the instructions given by tools/MK_GRADS/00README.txt. The MXE
package (Section 21.6) also has tools for visualization in the directory, postp/.

21.3.2 Restart for the SOM advection scheme
When the SOM tracer advection scheme is used (SOMADVEC), the model may read and write additional 9 restart files of the
2nd order moments for each tracer. Their names are BASENAME_SS_NN.YYYYMMDDHHMMSS, where BASENAME indicates acharacter string specified separately for input and out files by file_som_in and file_som_out in namelist nml_somadv,
SS elements of the moments (9 types), NN a tracer number (01 for temperature, 02 for salinity, numbers after 03 for othertracers), and the suﬃx for the date and time. The format of each file is the same as the tracer restart files. See Section 8.6
for the SOM scheme, and docs/README.Namelist for other items of nml_somadv.

21.3.3 Restart for the sea ice model
To handle restart files of the sea ice model, dedicated namelist blocks must be specified separately for input and output asfollows:An example namelist for sea ice restart input&nml_seaice_rst_in
file_icecat_restart_in = ’result/rs_ice’,/

An example namelist for sea ice restart output

&nml_seaice_rst_out
file_icecat_restart_out_temp = ’result/rs_ice’,nwrt_rst_ic= 480, !- set same as nstep_seaice/where file_icecat_restart_in and file_icecat_restart_out_temp specify the basenames of input and output
files, and nwrt_rst_ic specifies the time step interval to write instantaneous states (usually, set it the same as the time
step number of the integration). See docs/README.Namelist for other items. (Before MRI.COM ver.4.2, one must set
l_file_name_calendar = .true. in nml_seaice_time_conf to add a suﬃx of the date and time to the file name inthe same manner as the ocean model. After ver.4.4, the suﬃx is always added. )The restart file can be read by the following program.

– 266 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 21

User’s GuideProgram to read restart for the sea ice model

integer(4), parameter :: ncat = 5 ! number of thickness categories! ice concentration
real(8) :: aicen (1:imut,1:jmut,0:ncat), a0iceo(imut,jmut)! ice thickness
real(8) :: hicen (1:imut,1:jmut,0:ncat), hiceo (imut,jmut)! averaged sea-ice thickness
real(8) :: hin(1:imut,1:jmut,0:ncat), hi(imut,jmut)! snow depthreal(8) :: hsnwn (1:imut,1:jmut,0:ncat), hsnwo (imut,jmut)! averaged snow thickness
real(8) :: hsn(1:imut,1:jmut,0:ncat), hsnw (imut,jmut)! ice surface temperaturereal(8) :: tsfcin(1:imut,1:jmut,0:ncat), tsfci (imut,jmut)! ice temperaturereal(8) :: t1icen(1:imut,1:jmut,0:ncat)! sea surface skin temperaturereal(8) :: t0icen(1:imut,1:jmut,0:ncat)! sea surface skin salinity
real(8) :: s0n(1:imut,1:jmut,0:ncat)! skin temperature beneath the sea ice
real(8) :: t0iceo(imut,jmut)! skin temperature in the open leadsreal(8) :: t0icel(imut,jmut)! stress tensor
real(8) :: sigma1(imut,jmut), sigma2(imut,jmut), sigma3(imut,jmut)
integer(4),parameter :: nu = 10 ! device number
integer(4):: i = 0, m
open(nu, file=’result/rs_ice’, form=’unformatted’, && access=’direct’, recl=imut*jmut*8)
do m = 0, ncat ; i = i + 1 ; read(nu,rec=i) aicen(1:imut,1:jmut,m) ; enddo
do m = 0, ncat ; i = i + 1 ; read(nu,rec=i) hin (1:imut,1:jmut,m) ; enddo
do m = 0, ncat ; i = i + 1 ; read(nu,rec=i) hsn (1:imut,1:jmut,m) ; enddo
do m = 0, ncat ; i = i + 1 ; read(nu,rec=i) hicen(1:imut,1:jmut,m) ; enddo
do m = 0, ncat ; i = i + 1 ; read(nu,rec=i) hsnwn(1:imut,1:jmut,m) ; enddo
do m = 0, ncat ; i = i + 1 ; read(nu,rec=i) tsfcin(1:imut,1:jmut,m) ; enddo
do m = 0, ncat ; i = i + 1 ; read(nu,rec=i) t1icen(1:imut,1:jmut,m) ; enddo
do m = 0, ncat ; i = i + 1 ; read(nu,rec=i) t0icen(1:imut,1:jmut,m) ; enddo
do m = 0, ncat ; i = i + 1 ; read(nu,rec=i) s0n(1:imut,1:jmut,m) ; end do
i = i + 1 ; read(nu,rec=i) a0iceo(1:imut,1:jmut)
i = i + 1 ; read(nu,rec=i) hi(1:imut,1:jmut)
i = i + 1 ; read(nu,rec=i) hsnw(1:imut,1:jmut)
i = i + 1 ; read(nu,rec=i) hiceo(1:imut,1:jmut)
i = i + 1 ; read(nu,rec=i) hsnwo(1:imut,1:jmut)
i = i + 1 ; read(nu,rec=i) uice(1:imut,1:jmut)
i = i + 1 ; read(nu,rec=i) vice(1:imut,1:jmut)
i = i + 1 ; read(nu,rec=i) sigma1(1:imut,1:jmut)
i = i + 1 ; read(nu,rec=i) sigma2(1:imut,1:jmut)
i = i + 1 ; read(nu,rec=i) sigma3(1:imut,1:jmut)
close(nu)– 267 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201721.4 Execution

21.3.4 Control of input and output at run timea. Input

Input of restart files at run time is controlled by namelist listed on Table 21.5. By setting l_rst_in = .true. in
nml_run_ini_state, all necessary elements are read from restart files. However, the specification by l_rst_in may beoverridden for some elements by specifying namelists dedicated to those elements. Note that model tries to read restart
files for all passive tracers, thus namelist nmlrs_ptrc must be always specified appropriately.
Table21.5 Namelist blocks that control input of restart files.

namelist block
nml_run_ini_state

nml_barotropic_run
nml_baroclinic_run
nml_tracer_run
nml_vmix_run
nml_tide_run
nml_melyam_run
nml_nohkim_run
nml_gls_run
nml_density_run
nml_seaice_run
nml_tracer_data
nml_tracer_data

variable namel_rst_in

Usage.true. : read restart files and resume integration.false. (default): initial state is set internallyNo motion (u,v,ssh=0).
Initial condition for temperature and salinity are
deteremined by nml_tracer_run
l_rst_barotropic_in
default = l_rst_in
l_rst_baroclinic_in
default = l_rst_in
l_rst_tracer_in
default = l_rst_in
applicable to only temperature and salinity
l_rst_vmix_in
default = l_rst_in
l_rst_tide_in
default = l_rst_in
l_rst_melyam_in
default = l_rst_in
l_rst_nohkim_in
default = l_rst_in
l_rst_gls_in
default = l_rst_in
l_rst_density_in
default = l_rst_in
l_rst_seaice_in
default = .false.
adv_scheme%lrstin_som = .true. to read restart of SOM scheme
adv_scheme%lrstout_som = .true. to write restart of SOM schemeb. Output

Model always try to write restart files of the oceanic part according to the specification by the namelists listed on Table
21.4. Users have to specify those namelists appropriately to obtain restart files and terminate the model normally. Output
from the SOM scheme may be suppressed if adv_scheme%lrstout_som = .false. in nml_tracer_data.

21.4 Execution
To run a model, a shell script that handles input/output files, executes the compiled binary ogcm, and post-processes isusually prepared. The namelist parameters that control the time integration and the namelist file that specify sampling
have to be given.
• NAMELIST.name_model controls the job and gives parameters to sub-packages. See docs/README.Namelist for
details.
• NAMELIST.name_model.MONITOR specifies sampling. See docs/README.Monitor for details.
Parameters used by sub-packages are explained in corresponding chapters. Runtime parameters relevant to timeintegration are listed on Tables 21.6 through 21.11.

– 268 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 21

User’s Guide
Table21.6 Namelist nml_time_step.

variable namedt_secunits
secnstep_matsuno_interval

integeralpha_bryan_1984

real(8)l_monitor_time

logicaldescription
unit time interval for equation of motion
and tracer
time step interval with which Matsuno(Euler-backward) time-integration scheme
is usedacceleration parameter α of Bryan (1984).
See Section 2.3.1.
time step monitor is written to standard output ()usage
requireddefault = 12 steps
default = 1
default = .true.

Table21.7 Namelist nml_run_period.

variable namenstep_total

units
integer

description
total number of time step of this run

usage
required

Table21.8 Namelist nml_exp_start that specifies start time of the whole experiment.

variable name
year
month
day
hour
minutesecondunits
integer(4)integer(4)integer(4)integer(4)integer(4)integer(4)description

usage
default = −999
default = 1
default = 1
default = 0
default = 0
default = 0

year
month
day
hour
minutesecondTable21.9 Namelist nml_run_ini that specifies start time of this run.

variable name
year
month
day
hour
minutesecondunits
integer(4)integer(4)integer(4)integer(4)integer(4)integer(4)description

usage
default = 1
default = 1
default = 1
default = 0
default = 0
default = 0

year
month
day
hour
minutesecondTable21.10 Namelist nml_calendar that specifies treatment of leap year.

variable namel_force_leap

unitslogicall_leaplogicaldescription
Use l_leap to decide whether the current year
is leap or not.
the current year is a leap year or not.

usage
default = .false., the realistic calendar isfollowed
default = .false.

Table21.11 Namelist nml_stdout that specifies standard output (program log).

variable namel_stdout2file
file_base_stdout

units
logicall_debug

logicaldescription
Standard output is written to separate files for MPI processes
file name is name_model-file_base_stdout-stdout.XXXX(XXXX: mpi process number)debug mode

– 269 –

usage
default = .false.
default = .false.

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201721.5 Monitoring

21.5 Monitoring
This section summarizes outputs of the mean state (history) data used for monitoring and analyses.

21.5.1 History of the ocean and ice models
In MRI.COM version 4, specification of history outputs has been significantly changed and commonized in the oceanand sea ice models. Users prepare a namelist file named NAMELIST.name_model.MONITOR dedicated to history outputs(name_model can be specified by NAME_MODEL in configure.in, default=OGCM), and write the following namelist blocks,
nml_history, in it as needed,An example namelist for mean temperature output&nml_history
name= ’temperature’
file_base= ’result/hs_t’
interval_step = 48
[suffix= ’day’]/where name specifies the variable name to be output、file_base the basename of the history files, interval_step theoutput interval in terms of the integration time step. An option item of suffix specifies the depth of calendar date and timeused in the file suﬃx. In the above example, the time-mean temperature is written to the file, result/hs_t.YYYYMMDD, atevery 48-th time step.
History files are saved in a Fortran direct-access format, which can be read by a following program.
Program to read ocean model history datacharacter(8):: date = ’20010101’ !- for 1JAN2001
real(4):: d(imut,jmut,km)
integer(4),parameter :: nu = 10 ! device number
open(nu, file=’result/hs_t.’//date, form=’unformatted’, && access=’direct’, recl=imut*jmut*km*4)
read(nu,rec=1) dclose(nu)
A GrADS control file to visualize data is also made by default (result/hs_t.ctl in the above example).
The state variables that can be monitored by namelist nml_history are listed in docs/README.Monitor. There aremany options in nml_history, such as netCDF output, snapshot output, averaging in the model region, output in a
specified sub region, addition of an oﬀset value, multiplication by a factor, and so on. See docs/README.Monitor for theavailable options.
History outputs of the sea ice model are also specified by nml_history in the same manner as the ocean model. Statevariables that can be monitored are listed in the sea ice section of docs/README.Monitor. The format of output files
is also common, except that a missing value for grids of ocean without sea ice should be diﬀerent from the one for land
grids. By default, the value of 0.e0 is used for the former, while -9.99e33 for the latter. (These values can be changed by
nml_seaice_hst written in NAMELIST.name_model.)21.6 MRI.COM eXecution Environment (MXE)
We have been developing a package "MRI.COM eXecution Environment (MXE)" that aggregates programs for preprocessing, execution, postprocessing, and analysis of MRI.COM experiments. By using prepared tools, users can relatively
easily perform complex work of model building and various analyses. In fact, this package is also used as a common basis
for developing various ocean models at the Meteorological Research Institute.
The MXE package includes the following tools.
• Preprocessing tools for creating experiment setup files (directory prep)• Directory template for running MRI.COM (exp)• Postprocessing tools for visualizing output files (postp)• Various analysis tools (anl)– 270 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 21

User’s Guide• A Fortran library that provides basic subroutines such as grid and topography information (lib)The main programs are written in Fortran and sample shell scripts are also included for execution. Several tools in the
package have been confirmed by unit testing. This package is managed independently from MRI.COM, but it is provided
to users as an accompanying material. Since it is often updated like MRI.COM, it is recommended to use the latest version.
For details on how to use MXE, see the file README-MXE.md in the top directory (in Japanese).

21.7 Appendix
21.7.1 Model options
The model options are listed on Table 21.12. Only major options are listed here. Description about those related tobio-geochemical models can be found in Chapter 11. The description of all options for the latest version can be foundin src/README.Options. Though an expression like OGCM_PARALLEL is used in the source program, OGCM_ is omitted
when users specify options in configure.in.
Table21.12 Description of Model Options

Model option
BBL
BIHARMONICBULKKONDO2BULKNCAR

BULKITERCALALBSICALPPCARBONCBNHSTRUNCHLMA94CYCLICDIFAJSF2003
FSVISCGMANISOTROPGMTAPER

Description
uses the bottom boundary layer model
uses biharmonic operator for both horizontal viscosity and diﬀusion(*) If ISOPYCNAL is also selected, the biharmonic form is used only for viscosity and not for diﬀusion.
Kondo (1975) is used for the surface flux bulk formula.
Large and Yeager (2004; 2009) is used for the surface flux bulk formula. This option corresponds to
the COREs.
(*) BULKKONDO2, BULKNCAR is available only for HFLUX case.
Bulk transfer coeﬃcient is calculated using iterative method if the observed wind speed is not at
10 m.
(*) use with BULKKONDO2 and BULKNCAR
Sea-ice albedo is calculated using sea-ice conditions according to Los-Alamos model instead of using
a constant value
considers the time variation of pressure for the equation of state
bio-geochemical process is included(*) numtrc_p= 4 for Obata-Kitamura model; numtrc_p= 8 for NPZD model
atmospheric pCO2 is given from file(*) use with CARBON
shortwave penetration scheme with chlorophyll concentration by Morel and Antoine (1994)(*) use with NPZD and SOLARANGLE
uses zonally cyclic condition
sets large vertical diﬀusion coeﬃcient (1.0 m2 s−1 ) between unstable points instead of convective
adjustment
Program uses Fortran 2003 features
calculates viscosity explicitly in the barotropic momentum equation
Anisotropic horizontal variation of thickness diﬀusion is used(*) use with ISOPYCNALTaper GM vector stream function near the sea surface(*) cannot be used with SLIMIT, GMANISOTROP, AFC
Continued on next page– 271 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201721.7 Appendix

Model option
GMVAR
GLS
HFLUX
ICE
ICECLIMISOPYCNALISOTAPERLWDOWNMELYAM
MPDATAADVECNOHKIM
NOMATSUNONPZDOFFNESTPAR
OFFNESTSUB
PARALLEL
PARENT
RUNOFFSFLUXRSFLUXWSIDYN
SLIMITSLPSMAGHDSMAGORSOLARANGLE
SOMADVECSPHERICALSUBTAUBULK
TDEW
TIDE
TRCBIHARM

Table 21.12 – continued from previous page
Description
Horizontal thickness diﬀusion is allowed to vary in horizontal
Generic length scale model of Umlauf and Burchard (2003)
calculates sea surface heat flux using bulk formulasea ice is included
reading climatological sea-ice fractional area from fileuses isopycnal diﬀusion and Gent-McWilliams’ parameterization for eddy
induced tracer transport velocity (thickness diﬀusion)
Taper isopycnal diﬀusion coeﬃcient (ISOPYCNAL)
Longwave radiation data include only downward component instead of default net radiation
uses Mellor and Yamada Level 2.5 for mixed layer model
uses MPDATA for tracer advectionuses Noh’s mixed layer model
uses forward finite diﬀerence instead of Matsuno scheme
NPZD process is included(*) use with CARBON
Used as the lower-resolution part of an oﬀ-line one-way nesting calculation
Used as the higher-resolution part of an oﬀ-line one-way nesting calculation
parallel calculation using MPI. The number of zonally and meridionally partitioned regions shouldbe specified as NPARTX and NPARTY, respectively.
executed as low resolution model of the nesting calculation
uses river runoﬀ data(*) available only for WFLUX
SSS is restored to the climatological sea surface salinity as the salinity flux
calculates salinity flux converting from the surface freshwater flux(*) available only for WFLUX
sea-ice dynamics model with EVP rheology(*) available only for ICE
Tapers thickness diﬀusion near the sea surface
Sea surface is elevated/depressed according to surface atmospheric pressure
uses the Smagorinsky viscosity coeﬃcient multiplied by a constant ratio as the horizontal diﬀusion
coeﬃcient(*) available only for SMAGOR
uses the Smagorinsky parameterization for horizontal viscosity
solar insolation angle is considered in calculating shortwave penetration
uses second order moment advection by Prather (1986)
calculates scale factor semi-analytically for spherical coordinates
executed as a high resolution model of the nesting calculation
calculates the wind stress using bulk formulae by reading wind speed over the ocean
reads dew-point temperature and converts to specific humidity(*) available only for HFLUX
Tide producing forcing is activateduses biharmonic operator for horizontal diﬀusion(*) Should not be used with ISOPYCNAL
Continued on next page– 272 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 21

User’s Guide

Model option
TRIPOLARUTZQADVECVISANISO
VIS9P
VISBIHARM
VMBG3D
VVDIMP
WADJ
WFLUX
ZSTAR
MOVE
SCUP
SCUPCGCM
SCUPNEST

Table 21.12 – continued from previous page
Description
Tripolar system is used to construct model grids of a global model(*) Cannot be used with SPHERICAL"UTOPIA" and "QUICKEST" scheme can be used for horizontal and vertical tracer advection withultimate limiter
Anisotropic viscosity coeﬃcients are used(*) use with VIS9P
calculates the viscosity using adjacent 9 grid points
uses biharmonic operator for both horizontal viscosity
reads 3-D vertical viscosity and diﬀusion coeﬃcients from a filecalculates the vertical diﬀusion/viscosity by implicit method(*) it is automatically loaded if any mixed layer model is used or ISOPYCNAL option is selectedadjusts sea surface freshwater flux every time step to keep its global sum to be zero(*) available only for WFLUX
uses the sea surface freshwater flux to force the model
zstar (z*) vertical coordinate(*) KSGM = KM
used as ocean module for data assimilation (MOVE) systemuse simple coupler (SCUP) libraryused as an ocean module for a coupled model using scup for communication
on-line nesting(*) use with SCUP

– 273 –