TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Part IV

Evolution of tracers

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 8

Tracer advection schemes
This chapter describes tracer advection schemes available in MRI.COM. The tracer advection schemes available in
MRI.COM are: the weighted upcurrent scheme (Section 8.2), the Quadratic Upstream Interpolation for Convective
Kinematics (QUICK; Leonard, 1979) as described in Section 8.3 (QUICKADVEC option), a combination of the QUICK with
Estimated Streaming Terms (QUICKEST; Leonard, 1979) for vertical advection (Section 8.4) and the Uniformly ThirdOrder Polynomial Interpolation Algorithm (UTOPIA; Leonard et al., 1993), which is a two-dimensional generalization of
the QUICKEST scheme, for horizontal advection (Section 8.5) with UTZQADVEC option. For QUICKEST and UTOPIA,
a flux limiter that prevents unrealistic extrema is applied (Leonard et al., 1994).
The above schemes seek to improve the accuracy by refining the finite-diﬀerence expression at the cell faces. There areother approaches. The Second Order Moment (SOM; Prather, 1986) scheme seeks to improve the accuracy by considering
the distribution within the cell and is available in MRI.COM through SOMADVEC option (Section 8.6). The Multidimensional
Positive Definite Advection Transport Algorithm (MPDATA; Smolarkiewicz and Margolin, 1998) scheme seeks to cancel
the numerical diﬀusion of the lower order scheme by using anti-diﬀusive velocities and is available in MRI.COM through
MPDATAADVEC option (Section 8.7).
Diﬀerent advection schemes can be used for individual tracers. They should be chosen from among compiled schemes
at run time (Chapter 13).

8.1 Finite volume or flux form method
The equation governing the time change of tracer T is{}∂(z s hψ uT ) ∂(z s hµ vT )∂(z s T )1∂(z s ṡT ) ∂(z s T )+++=+ A(T ) = −z s ∇ · FT + z s QT ,∂th µ hψ
∂µ
∂ψ
∂s
∂t

(8.1)

where A(T ) represents the advection operator for T. The finite diﬀerence form of the l.h.s. of (8.1) is given by firstconsidering a control cell volume and then calculating fluxes through each cell face, and setting their divergence andconvergence to be the time change rate at the grid cell (Figure 8.1). In finite diﬀerence form, this is expressed as follows:
Ti,n+1∆Vi,n+1= Ti,n−1∆Vi,n−1j,k− 1j,k− 1j,k− 1j,k− 122

2

2+ 2 ∆t {FXAi− 1 , j,k− 1 − FXAi+ 1 , j,k− 1 + FYAi, j− 1 ,k− 1 − FYAi, j+ 1 ,k− 1 + FZAi, j,k − FZAi, j,k−1 }, (8.2)
2

2

2

2

2

2

2

2

where ∆V is the volume of the grid cell, and FXA, FYA, and FZA represent (flux due to advection) × (area of the cellboundary). The same consideration may be applied for the discretization of the diﬀusion term.
As explained in Chapter 3, ∆V varies with time and is given by (3.28),∆Vi, j,k− 1 ≡ (volt)i, j,k− 1 = (volu_bl)i+ 1 , j+ 1 ,k− 1 + (volu_tl)i+ 1 , j− 1 ,k− 1
2

2

2

2

2

2

2

2

+ (volu_br)i− 1 , j+ 1 ,k− 1 + (volu_tr)i− 1 , j− 1 ,k− 1 .
2

2

2

2

2

2

(8.3)

The volume of the left-lower quarter of the T-cells at (i, j) (corresponding to the southwestern part in geographic
coordinates) is represented by (volt_tr)i− 1 , j− 1 . Similarly, (volt_tl)i+ 1 , j− 1 is the right-lower (southeastern),2222(volt_br)i− 1 , j+ 1 is the left-upper (northwestern), and (volt_bl)i+ 1 , j+ 1 is the right-upper (northeastern) quarter of a2222T-cell at (i, j).

– 85 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.2 Weighted upcurrent scheme
Fluxes due to advection are given as follows:
T
FXAi+ 1 , j,k− 1 = Ui+
T 1
1,1, j,k− 1 i+ , j,k−

(8.4)

FYAi, j+ 1 ,k− 1 = Vi,Tj+ 1 ,k− 1 Ti, j+ 1 ,k− 1 ,

(8.5)

2

2

2

2

2

2

2

2

2

2

2

2

FZAi, j,k = Wi,Tj,k Ti, j,k ,(8.6)

where horizontal volume transport U T and V T are defined as follows, using (5.13) to (5.18):
T
Ui+
=
1
, j,k− 12

2

Vi,Tj+ 1 ,k− 1 =22∆yi+ 1 , j (2

2)
ui+ 1 , j− 1 ,k− 1 ∆zi+ 1 , j− 1 ,k− 1 + ui+ 1 , j+ 1 ,k− 1 ∆zi+ 1 , j+ 1 ,k− 1 ,

(8.7)

)
vi− 1 , j+ 1 ,k− 1 ∆zi− 1 , j+ 1 ,k− 1 + vi+ 1 , j+ 1 ,k− 1 ∆zi+ 1 , j+ 1 ,k− 1 .

(8.8)

2

∆x i, j+ 1 (
2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

Vertical volume transport W T is then obtained by diagnostically solving (5.12). Moreover, the vertical velocity w, which
is necessary for using QUICKEST, is calculated as follows (w is not needed except for QUICKEST):
Wi,Tj,k = wi, j,k × (areat)i, j,k− 1 ,

(8.9)

2

where
(areat)i, j,k− 1 = (a_tr)i− 1 , j− 1 × (aexl)i− 1 , j− 1 ,k− 1 + (a_tl)i+ 1 , j− 1 × (aexl)i+ 1 , j− 1 ,k− 1
2

2

2

2

2

2

2

2

2

2

2

+ (a_br)i− 1 , j+ 1 × (aexl)i− 1 , j+ 1 ,k− 1 + (a_bl)i+ 1 , j+ 1 × (aexl)i+ 1 , j+ 1 ,k− 1 .
2

2

2

2

2

2

2

2

2

2

(8.10)An array aexl is set to be unity if the corresponding U-cell is a sea cell and set to be zero otherwise.
This formulation does not depend on the choice of the advection scheme. The diﬀerence arises from the way ofdetermining cell boundary values of tracer, Ti+ 1 , j,k− 1 , Ti, j+ 1 ,k− 1 , and Ti, j,k . There are several choices for tracer advection22
22
schemes, which are explained in the following sections.
• Weighted upcurrent scheme (adv_scheme%name = ”upc”). This scheme is always available for use. Specify theweighting ratio listed on Table 8.1.
• Quadratic Upstream Interpolation for Convective Kinematics (QUICK; Leonard, 1979) for both horizontal andvertical direction (QUICKADVEC option; adv_scheme%name = ”quick”).
• The combination of the Uniformly Third-Order Polynomial Interpolation Algorithm (UTOPIA; Leonard et al.,1993) for horizontal direction and QUICK with Estimated Streaming Terms (QUICKEST; Leonard, 1979) for thevertical direction (UTZQADVEC option; adv_scheme%name = ”utzqadv”).
• Second Order Moment (SOM; Prather, 1986) schemes for both horizontal and vertical direction (SOMADVEC option;
adv_scheme%name = ”som”).
• Multidimensional Positive Definite Advection Transport Algorithm (MPDATA; Smolarkiewicz and Margolin, 1998)
for all directions (MPDATAADVEC option; adv_scheme%name = ”mpdata”).

8.2 Weighted upcurrent scheme
In the weighted upcurrent scheme (e.g., Suginohara and Aoki, 1991; Yamanaka et al., 2000), the cell boundary value isdetermined by a weighted average of the upcurrent scheme and the centered finite diﬀerence scheme.
The upcurrent scheme employs the upstream value as the cell boundary value: Ti, j,k− 12 ,upcurrent= Ti+1, j,k− 1 ,i+ 12 , j,k− 12

T



2

if
if

Ui+ 1 , j,k− 1 > 0,2
2
Ui+ 1 , j,k− 1 < 0.
2

(8.11)

2

The centered finite diﬀerence scheme uses the average between the two neighboring points of tracer as the cell boundary
value:
Ti+1, j,k− 1 + Ti, j,k− 122
center.
(8.12)Ti+
1
1 =2
2 , j,k− 2– 86 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 8

Tracer advection schemes

V Ti, j+1/2
UVi+1/2, j+1/2

UVi-1/2, j+1/2
U Ti-1/2, j

U Ti+1/2, jTi, jUVi-1/2, j-1/2

UVi+1/2, j-1/2
V Ti, j-1/2
W Ti, kUVi-1/2, k+1/2

Ti, k+1/2

UVi+1/2, k+1/2

W Ti, k+1
Figure8.1 Grid arrangement around TS-Box (Upper: Views from the upper, Below: views from the horizontal).
Fluxes represented by an arrow are calculated.
Table8.1 Namelist nml_tracer_adv for weighted upcurrent scheme

variable nameweight_upcurrent_horz

units
1

weight_upcurrent_vert

1

description
upstream-side weighting ratio for thehorizontal advection
upstream-side weighting ratio for thevertical advectionusage
1.0: up-current, 0.5: centered diﬀerence
1.0: up-current, 0.5: centered diﬀerence

Taking the ratio of the upcurrent scheme to be α (0 ≤ α ≤ 1), the tracer flux at the eastern face of a grid cell is expressed
as follows:
UT 1+ UT 1UT 1− UT 1
i+ 2 , j,k− 12
i+ 2 , j,k− 12
i+ 2 , j,k− 21
i+ 2 , j,k− 12T+T
FXAi+ 1 , j,k− 1 =α 
i
i+1


2
2
2
2




Ti+1, j,k− 1 + Ti, j,k− 1
22T+ (1 − α)Ui+
1
1
22 , j,k− 2
UT 1
UT 1



1(
(


i+ 2 , j,k− 12 )i+ 2 , j,k− 21 )1T.
= Ui+ 1 , j,k− 1 
1
+
α
T
+
1
−
α
T
i
i+1
T
T


2

2
2 2
U 1
U 1
1
1
i+
,
j,k−
i+
,
j,k−


2
2
2
2

(8.13)(8.14)

Mainly for the sake of computational eﬃciency, we give a parameter β = 12 (α + 1) instead of α at run time. Diﬀerent
parameters for horizontal and vertical directions may be given, which are listed on Table 8.1.

– 87 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.3 QUICK scheme

8.3 QUICK scheme
This section explains the QUICK scheme. In the QUICK scheme, the cell boundary value is interpolated by a quadratic
function, using three points, with one of them added from the upstream side (Figure 8.2).

T
Ti+1

quadratic function thatpasses three points

Ti

Ti+1/2

Ti-1
Ui+1/2

i-1

X

i+1

i
T-box

Figure8.2 Schematic for interpolation: Ti represents the value of a tracer at T-point with index i, and Ti+ 1 is the cell2boundary value. In the QUICK scheme, Ti+ 1 is interpolated by a quadratic function that passes Ti and the neighboring
2
T-point values, Ti−1 and Ti+1 .

Originally, cell boundary values in the QUICK scheme are given as follows:∆x i Ti+1, j,k− 1 + ∆x i+1Ti, j,k− 1∆x i+1 ∆x i
ci+ 1 , j,k− 1 ,2
2
4∆y j+1 ∆y j
2
2
Ti, j+ 1 ,k− 1 =−d i, j+ 1 ,k− 1 ,22
22∆y j+1 + ∆y j4∆zk− 1 Ti, j+1,k+ 1 + ∆zk+ 1 Ti, j,k− 1 ∆zk− 1 ∆zk+ 1
2
2
2
2
2
2
Ti, j,k =−ei, j,k ,∆zk+ 1 + ∆zk− 1
4
Ti+ 1 , j,k− 1 =
2

2

2

2

∆x i+1 + ∆x i
∆y j Ti, j+1,k− 1 + ∆y j+1Ti, j,k− 12−(8.15)(8.16)(8.17)

2

where c, d, and e are defined depending on the direction of the mass flux as follows:ci+ 1 , j,k− 1
2

2

=
=

d i, j+ 1 ,k− 1
2

2

=
=

ei, j,k==∆x i δ x δ x Ti, j,k− 12(≡ cp ),x2∆x i∆x i+1 δ x δ x Ti+1, j,k− 1
2(≡ cm ),x2∆x i+1∆y j δ y δ y Ti, j,k− 12(≡ d p ),y2∆y j∆y j+1 δ y δ y Ti, j+1,k− 12(≡ d m ),y2∆y j+1∆zk+ 1 δz δz Ti, j,k+ 1
2
2(≡ e p ),z
2∆zi, j,k+ 1
2∆zk− 1 δz δz Ti, j,k− 12
2(≡ em ),z2∆zi, j,k− 1
2– 88 –

if

UT 1

> 0,

if

UT 1

< 0,

if

VT

> 0,

if

VT

< 0,

if

Wi,Tj,k > 0,ifWi,Tj,k < 0.

i+ 2 , j,k− 12
i+ 2 , j,k− 12

i, j+ 12 ,k− 21

i, j+ 12 ,k− 21(8.18)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 8

Tracer advection schemes

The finite diﬀerence operators are defined as follows (definitions in y and z directions are the same):δ x Ai ≡x

Ai ≡Ai+ 1 − Ai− 1
2

2

∆x i
Ai+ 1 + Ai− 1
2

2

2

,

Ai+1 − Ai,∆x i+ 1

δ x Ai+ 1 ≡
2

2

,

x

Ai+ 1 ≡
2

Ai+1 + Ai.
2(8.19)

Letting cp , d p , and e p represent their values for positive velocity at the cell boundary and cm , d m , and em represent theirvalues for negative velocity at cell boundary and taking
ca = cm + cp
cd = cm − cp
da = dm + d p
dd = dm − d p
e a = em + e p
e d = em − e p,

(8.20)(8.21)(8.22)(8.23)(8.24)(8.25) ∆x i, j T∆x i+1, j ∆x i, j
i+1, j,k− 21 + ∆x i+1, j Ti, j,k− 12T
FXAi+ 1 , j,k− 1 = Ui+
c
−
1
1
a
1
1 
i+ 2 , j,k− 222∆x i+1, j + ∆x i, j
8
2 , j,k− 2 ∆x i+1, j ∆x i, j
T
+ |Ui+
cd i+ 1 , j,k− 1 ,
1|
1
2
2
8
2 , j,k− 2 ∆yi, j T
1 + ∆yi, j+1T
∆yi, j+1 ∆yi, j
i, j+1,k− 2i, j,k− 12
d a i, j+ 1 ,k− 1 −FYAi, j+ 1 ,k− 1 = Vi,Tj+ 1 ,k− 1 2
2
22∆yi, j+1 + ∆yi, j
8
2
2 



(8.26)

we obtain∆yi, j+1 ∆yi, j
d d i, j+ 1 ,k− 1 ,
2
2
8
2
2
 ∆z

∆zi, j,k+ 1 ∆zi, j,k− 1
i, j,k− 21 Ti, j,k+ 12 + ∆zi, j,k+ 12 Ti, j,k− 21
2
2
FZAi, j,k = Wi,Tj,k eai, j,k −∆zi, j,k+ 1 + ∆zi, j,k− 1
8


2
2
∆zi, j,k+ 1 ∆zi, j,k− 122
ed i, j,k .
+ |Wi,Tj,k− 1 |82+ |Vi,Tj+ 1 ,k− 1 |(8.27)(8.28)

Equation (8.26) can be rewritten as
T
FXAi+ 1 , j,k− 1 ≃ Ui+
1 T̃i+ 1 , j,k− 1 + AQ
1
2
2
2
2
2 , j,k− 2∂ 3Ti+ 1 , j,k− 1
2

∂ x3

2

,

(8.29)

where T̃i+ 1 , j,k− 1 is the value of T at the cell boundary interpolated by the cubic polynomial, and22T
AQ = |Ui+|
1, j,k− 12∆x i+1 ∆x i+ 1 ∆x i

2

2

8

.

(8.30)

Although the time integration for advection is done by the leap-frog scheme, the second term on the r.h.s. of (8.29) has
a biharmonic diﬀusion form, and thus the forward scheme is used to achieve calculation stability (Holland et al., 1998).
A similar procedure is applied for the north-south and vertical directions.
The weighted up-current scheme is used for vertical direction if wi, j,k−1 > 0 and the T-point at (i, j, k + 12 ) is below
the bottom. The upstream-side weighting ratio is given by the user as the namelist parameter specified for the up-current
scheme (Table 8.1).

8.4 QUICKEST for vertical advection
This section describes the specific expression and the accuracy of the QUICK with Estimated Streaming Terms (QUICKEST; Leonard, 1979) for vertical advection.

– 89 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.4 QUICKEST for vertical advection
Consider a one-dimensional equation of advection for incompressible fluid∂T∂+(wT ) = 0,∂t∂z(8.31)

where w is a constant. Although the velocities are not uniform in the real three dimensional ocean, we assume a constantvelocity for simplicity.
Following the notation of vertical grid points and their indices (Section 3.2), tracers are defined at the center (k − 12 )of the vertical cells and vertical velocities are defined at the top (k − 1) and bottom (k) faces of the vertical cells. The
following relation holds for the vertical grid spacing:∆zk =∆zk+ 1 + ∆zk− 1
2

2

2

.

(8.32)In QUICKEST, the distribution of tracer T is defined using the second order interpolations, and the mean value during
a time step at the cell face (boundary of two adjacent tracer cells) is calculated. The coeﬃcients for the second orderinterpolation are calculated first. A Taylor expansion of T about point zk gives
2∆zk− 3∆zk− 32
2+ ∆zk− 1 + + c2 *+ ∆zk− 1 + + O(∆z 3 ),Tk− 3 = c0 + c1 *
2
2
2
, 2
, 2
2
∆z
∆zk− 1k− 21
2+ c2+ O(∆z 3 ),Tk− 1 = c0 + c1
2
2
4
∆z 2 1∆zk+ 1k+ 2
2+ c2+ O(∆z 3 ),Tk+ 1 = c0 − c122
42
∆zk+ 3∆zk+ 32
2
Tk+ 3 = c0 − c1 *+ ∆zk+ 1 + + c2 *+ ∆zk+ 1 + + O(∆z 3 ).
2
2
2
, 2
, 2
-

(8.33)(8.34)(8.35)(8.36)

Coeﬃcients c0 , c1 , and c2 can be solved using three of the four equations (8.33), (8.34), (8.35), and (8.36). The three
upstream-side equations are chosen. When w > 0 (w < 0), equations (8.34), (8.35), and (8.36) ((8.33), (8.34), and (8.35))
are used. The solution is as follows:
c0 =
c1 =

Tk− 1 ∆zk+ 1 + Tk+ 1 ∆zk− 1
2

2

2

2∆zk

Tk− 1 − Tk+ 1
2

2

∆zk

−
(T

2

−

∆zk+ 1 ∆zk− 12∆zk− 1 − ∆zk+ 1
2

2

4
2

2

c2,

c2,

(8.37)(8.38))−Tk+ 1
Tk+ 1 −Tk+ 3
k− 1

1
2
2
2
2

−


∆z k +∆z k+1∆z k∆z k+1()
c2 = 
Tk− 3 −Tk− 1
Tk− 1 −Tk+ 1

1
2
2
2
2

 ∆z +∆z−∆z k−1∆z k k−1 k(w > 0),(8.39)(w < 0).

Next, equation (8.31) is integrated over one time step and one grid cell.
∫ t n+1∫ zk−1
dt

tn

dz
zk

The r.h.s. of (8.40) can be written as
−

∂T
=−
∂t

∫ t n+1
tn∫ t n+1∫ zk−1
dt

tn

dz
zk

∂
(wT ).
∂zdt(wu Tu − wl Tl ),(8.40)(8.41)

where subscript u (l) denotes z = zk−1 (z = zk ). Assuming that w does not depend on time,∫ t n+1
tn

dtTl =∫ 0−wl ∆t[c0n + c1n ξ + c2n ξ 2 + O(∆z 3 )]

dξ
.
wl

(8.42)

Thus expression (8.41) becomes−∆t(wu Tun − wl Tln ) + O(∆z 3 w∆t),– 90 –(8.43)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 8

Tracer advection schemes

where
1
wl ∆t

∫ 0

(c0n + c1n ξ + c2n ξ 2 )dξ−wl ∆tcn
cn= c0n − 1 wl ∆t + 2 wl2 ∆t 2 .

Tln =

2

3

(8.44)

Using up to the second order terms of a Taylor expansion, the l.h.s. of (8.40) can be written as follows:∫ t n+1∫ zk−1
dt

tn

zk



∆z 2 1
k− 2
∂T

n+1
n+1
n
n
3 
dz
= ∆zk− 1 Tk− 1 − Tk− 1 +(Tzz k− 1 − Tzz k− 1 ) + O(∆z )  ,
2
∂t
24
2
2
2
 2


(8.45)

where
n

∂Tzz
+ O(∆t 2 )∂t k− 1
2
2
]n
2 [
∂
∂
+ O(∆t 2 )= −∆t 2
(wT )
1
∂z ∂z
k−

n
n+1
= ∆t
Tzz k−
1 − Tzz
k− 1
2

2

∂
n
2
= −∆t (wTzz )k−1 + O(∆t )∂z2 wu Tzz un − wl Tzz ln  + O(w∆t∆z) + O(∆t 2 ).
= −∆t ∆zk− 1
2(8.46)The expression for the r.h.s. of (8.45) becomes∆z 2 1 w T n − w T n
k− 2
u zz ul zz ln
3 n+1∆zk− 1 Tk− 1 − Tk− 1 −∆t+ O(∆z )  + O(w∆t∆z 4 ) + O(∆z 3 ∆t 2 ).
224∆zk− 12
 2

2
Based on (8.43) and (8.47), the discretized forecasting equation is expressed as follows:
∆t
n+1
n
Tk−
−
1 = T
k− 12
∆zk− 1
2
2



∆z 2
w T n − w T n − k− 21 (w T n − w T n )  + O(α∆z 4 ) + O(∆z 2 ∆t 2 ),l l
u zz ul zz l  u u24


(8.47)(8.48)

where
w∆t
< 1,
∆z
n
Tzz l = 2c2 + O(∆z).
α≡(8.49)(8.50)

The accuracy of equation (8.48) is max(O(∆z 4 ), O(∆z 2 ∆t 2 )).

8.5 UTOPIA for horizontal advection
The Uniformly Third Order Polynomial Interpolation Algorithm (UTOPIA; Leonard et al., 1993) is an advection schemethat can be regarded as a multi-dimensional version of QUICKEST. In UTZQADVEC option, horizontally two-dimensional
advection is calculated using UTOPIA. Vertical advection is calculated separately using QUICKEST.
Since grid intervals could be variable in both zonal and meridional directions in MRI.COM, UTOPIA is formulated
based on a variable grid interval. It is assumed that the tracer cell is subdivided by the borderlines of the velocity cells
into four boxes with (almost) identical area.
Consider an equation of advection:
1 ∂
1 ∂
∂T
+
(hψ uT ) +(hµ vT ) = 0.
∂th µ hψ ∂ µhµ hψ ∂ψ
Integrated over a tracer cell and for one time step,∫ ψL +∆ψL /2∫ µ L +∆µ L /2dψdµ( χ n+1 − χ n ) = −∆t(urnTrn ∆yr − ulnTln ∆yl + vunTun ∆x u − vdnTdn ∆x d ),
ψL −∆ψL /2µ L −∆µ L /2

– 91 –(8.51)(8.52)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.5 UTOPIA for horizontal advection
where χ ≡ hµ hψ T and Trn etc. on the r.h.s. are the face values described later. On the l.h.s. of (8.52), the second-orderinterpolation of χ is used to integrate the terms. The Taylor expansion of χ about L is given as follows (see Figure 8.3 forthe label of the point):
χ = χ L + a10 (µ − µ L ) + a20 (µ − µ L ) 2 + a01 (ψ − ψ L ) + a02 (ψ − ψ L ) 2 + a11 (µ − µ L )(ψ − ψ L ).

(8.53)

Then values at points E, W, N, and S are

NW

N

W

L

E

SW

S

SE

WW

u

r

SS
l

d

Figure8.3 Labels of tracer grid points (upper case characters) and faces (lower case characters).

χ E = χ L + a10 ∆µr + a20 ∆µ2r ,(8.54)

χ L − a10 ∆µl + a20 ∆µ2l ,χ L + a01 ∆ψu + a02 ∆ψu2 ,χ L − a01 ∆ψ d + a02 ∆ψ d2 ,(8.55)

χW =
χN =
χS =

(8.56)(8.57)

where∆ψ L + ∆ψ N,
2
∆ψ L + ∆ψS∆ψ d ≡,
2∆µ L + ∆µ E
∆µr ≡
,
2
∆µ L + ∆µW
.
∆µl ≡
2
∆ψu ≡

(8.58)(8.59)(8.60)(8.61)

Using these known values, the following parameters are obtained,χE − χLχ L − χW
+ ∆µr
∆µr
∆µl
a10 =
,
∆µr + ∆µl
χE − χLχ L − χW
−
∆µr
∆µl
a20 =
,
∆µr + ∆µl
χN − χLχ L − χS
∆ψ d
+ ∆ψu
∆ψu
∆ψ d
a01 =
,
∆ψu + ∆ψ d
χN − χLχ L − χS
−
∆ψu
∆ψ d
a02 =
.
∆ψu + ∆ψ d∆µl

– 92 –(8.62)(8.63)(8.64)(8.65)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 8

Tracer advection schemes

Substituting (8.53) into the l.h.s. of (8.52) yields
∆ψ 2L n+1∆µ2L n+1
n
n 
n
)
+
)
−
a
−
a
∆µ L ∆ψ L  χ n+1
−
χ
+
(a
(a
L
20
02  .
12 20
12 02

 L

(8.66)

Using equation (8.51), the following approximation is allowed:nn  hψr unT n − hψl unT n
hµu vunTµµu− hµd vdnTµµd

r µµr
l µµl
n
n+1
a20
− a20
= −∆t 
+
 ,
∆µ L
∆ψ L


n n
n n
 hψr unT n − hψl unT n

T
h
v
T
−
h
v
µu
µd
r ψψr
u ψψu
l ψψl
d ψψd 
n+1
n
a02
− a02
= −∆t 
+
 ,
∆µ L
∆ψ L



(8.67)(8.68)

n is the value of the second order derivative at the right face r, whose expression is similar to that of c described
where Tµµr
20later. Therefore, under a suitable approximation,

TLn+1 = TLn −∆t n ˜n(u T ∆yr − ulnT˜ln ∆yl + vunT˜un ∆x u − vdnT˜dn ∆x d ),∆SL r r(8.69)

where∆µ2L n∆ψ 2L nTµµl −T ,24
24 ψψl∆ψ 2L n∆µ2L nTµµd −T.
T˜dn = Tdn −2424 ψψdT˜ln = Tln −(8.70)(8.71)

Next, the expressions for Tln and Tdn are required. The term Tln is the average over the hatched area of Figure 8.4, andthe values of T n are given as the second order interpolation about l of Figure 8.3. Similar operations will be used to obtain
the expression for Tdn .

W

L

vl∆t

ul∆t
Figure8.4 Area used to average tracer values for the face lFirst, Taylor expansions of T n about l and d are written as follows:
T n |l =c00 + c10 (µ − µl ) + c20 (µ − µl ) 2 + c01 (ψ − ψ L ) + c02 (ψ − ψ L ) 2 + c11 (µ − µl )(ψ − ψ L ),(8.72)

T |d =d 00 + d 10 (µ − µl ) + d 20 (µ − µl ) + d 01 (ψ − ψ L ) + d 02 (ψ − ψ L ) + d 11 (µ − µl )(ψ − ψ L ).

(8.73)n2

2– 93 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.5 UTOPIA for horizontal advection
The T n values at eight points around l are,∆µ2∆µ L
+ c20 L ,
2
4
2
∆µ
∆µ
W
n
+ c20 W ,TW
= c00 − c10
2
4
(
)
(
)2
∆µ E
∆µ E
n
TE = c00 + c10 ∆µ L ++ c20 ∆µ L +
,
2
2
(
)
(
)2
∆µWW
∆µWW
n
TWW
= c00 − c10 ∆µW ++ c20 ∆µW +
,
2
2
∆µ L
TNn = TLn + c01 ∆ψu + c02 ∆ψu2 + c11
∆ψu,
2
∆µ L
TSn = TLn − c01 ∆ψ d + c02 ∆ψ d2 − c11
∆ψ d,
2
∆µW
n
n
∆ψu,
TNW
= TW
+ c01 ∆ψu + c02 ∆ψu2 − c11
2
∆µW
n
n
TSW
= TW
− c01 ∆ψ d + c02 ∆ψ d2 + c11∆ψ d .
2
TLn = c00 + c10(8.74)(8.75)(8.76)(8.77)(8.78)(8.79)(8.80)(8.81)

To obtain all six coeﬃcients, six of these equations (points) are used. The equations are chosen according to the following
flow direction.
uln > 0,vln > 0 ⇒ L, W, WW, S, NW, SW(8.82)

uln < 0,uln > 0,uln < 0,vln > 0 ⇒ L, W, E, N, S, SW
vln < 0 ⇒ L, W, WW, N, NW, SW
vln < 0 ⇒ L, W, E, N, S, NW(8.83)(8.84)(8.85)

From equations (8.74) and (8.75),
n∆µW TLn + ∆µ L TW∆µ L ∆µW
− c20
,
2∆µl
4
n
T n − TW∆µ L − ∆µW
c10 = L
− c20
.
∆µl
2

c00 =

(8.86)(8.87)

When uln > 0, from (8.74) and (8.77),
nn
TLn − TW
T n − TWW
− W
∆µl
∆µll
c20 =
,
∆µl + ∆µll∆µW + ∆µWW.
where ∆µll ≡2

Using equations (8.80) and (8.81),

When uln < 0, from (8.75) and (8.76),

Using equations (8.78) and (8.79),nnn
T n − TSW
TNW
− TW
− W
∆ψu
∆ψ d
c02 =
.
∆ψu + ∆ψ dnTEn − TLn TLn − TW
−
∆µr
∆µl
c20 =
.
∆µr + ∆µl

TNn − TLn TLn − TSn
−
∆ψu
∆ψ d
.
c02 =
∆ψu + ∆ψ d
– 94 –(8.88)(8.89)(8.90)(8.91)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 8

Tracer advection schemes

When vln > 0, from (8.79) and (8.81),
nn∆µW (TLn − TSn ) + ∆µ L (TW− TSW)+ c02 ∆ψ d,
2∆µl ∆ψ dn
T n − TW
− TSn + TLn
c11 = SW.
∆µl ∆ψ dc01 =(8.92)(8.93)

When vln < 0, from (8.78) and (8.80),nn∆µW (TNn − TLn ) + ∆µ L (TNW− TW)
− c02 ∆ψu,2∆µl ∆ψunn
T n − TLn − TNW+ TW
c11 = N.
∆µl ∆ψuc01 =(8.94)(8.95)

Next, using equation (8.73), the T n values at eight points around d are
∆ψ 2
∆ψ L
+ d 02 L ,
2
4
∆ψ 2
∆ψS
TSn = d 00 − d 01
+ d 02 S ,
2
4
)
)2
(
(
∆ψ N
∆ψ N
+ d 02 ∆ψ L +,TNn = d 00 + d 01 ∆ψ L +
2
2
(
)
(
)2
∆ψSS
∆ψSS
n
TSS
= d 00 − d 01 ∆ψS ++ d 02 ∆ψS +
,
2
2
∆ψ L
TEn = TLn + d 10 ∆µr + d 20 ∆µ2r + d 11
∆µr ,
2
∆ψ L
n
TW
= TLn − d 10 ∆µl + d 20 ∆µ2l − d 11
∆µl,
2
∆ψS
n
= TSn + d 10 ∆µr + d 20 ∆µ2r − d 11
TSE
∆µr ,
2
∆ψS
n
= TSn − d 10 ∆µl + d 20 ∆µ2l + d 11
TSW
∆µl .
2
TLn = d 00 + d 01

(8.96)(8.97)(8.98)(8.99)(8.100)(8.101)(8.102)(8.103)

From equations (8.96) and (8.97),∆ψS TLn + ∆ψ L TSn∆ψ L ∆ψS− d 02,2∆ψ d
4
T n − TSn∆ψ L − ∆ψSd 01 = L.
− d 02∆ψ d
2d 00 =(8.104)(8.105)

When vdn > 0, from (8.96) and (8.99),nTLn − TSn TSn − TSS
−
∆ψ d
∆ψ dd
d 02 =,∆ψ d + ∆ψ dd∆ψS + ∆ψSS
where ∆ψ dd ≡.
2

From (8.102) and (8.103),

When vdn < 0, from (8.97) and (8.98),nn
TSE− TSn TSn − TSW
−
∆µr
∆µl
.
d 20 =∆µr + ∆µl

TNn − TL TLn − TSn
−
∆ψu
∆ψ d
d 02 =.
∆ψu + ∆ψ d
– 95 –(8.106)(8.107)(8.108)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.5 UTOPIA for horizontal advection
From (8.100) and (8.101),nTEn − TLn TLn − TW
−
∆µr
∆µl
d 20 =.
∆µr + ∆µl

(8.109)When udn > 0, from (8.101) and (8.103),nn∆ψS (TLn − TW
) + ∆ψ L (TSn − TSW)+ d 20 ∆µl,2∆ψ d ∆µlnn
T n − TSn − TW
+ TSW
d 11 = L.
∆ψ d ∆µl

d 10 =(8.110)(8.111)

When udn < 0, from (8.100) and (8.102),n∆ψS (TEn − TLn ) + ∆ψ L (TSE− TSn )
− d 20 ∆µr ,2∆ψ d ∆µrn
T n − TSE
− TLn + TSn.
d 11 = E∆ψ d ∆µr

d 10 =(8.112)(8.113)The value of Tln is the average of T n over the hatched area of Figure 8.4. Defining
ξln =

uln
,
hµl

vln
,
hψl

η ln =(8.114)

we have∫ ψL +∆ψL /2 ∫ µl
1nT n dµdψ
Tl = nξl ∆t∆ψ L  ψL −∆ψL /2 µl −ξln ∆t+∫ ψL −∆ψL /2

ψL −∆ψL /2−ηln ∆t∫ µl + ξlnn (ψ−ψL +∆ψL /2)η[lµl −ξln ∆t{T (ψ) − T (ψ + ∆ψ L )}dµdψ 

]
n

n

1
1
1
= c00 − η ln ∆t c01 +∆ψ 2 + ∗ (η n ∆t) 2 c02
212 L 3 l1
11− ξln ∆t c10 + (ξln ∆t) 2 c20 + † ξln η ln ∆t 2 c11 .
2
3
3(8.115)This is the result for uln > 0 and vln > 0. The result is the same independent of the sign of uln and vln .
Similarly,
11
Tdn = d 00 − η dn ∆t d 01 + (η dn ∆t) 2 d 02
2
[ 3
]
1
1
1
1 n
∆µ2L + ‡ (ξ dn ∆t) 2 d 20 + § ξ dn η dn ∆t 2 d 11,− ξ d ∆t d 10 +
2
12
3
3
where
ξ dn =udn,
hµd

η dn =∗ 1 in COCO2
† 1 in COCO4
‡ 1 in COCO
2
1§
4 in COCO– 96 –

vdn.
hψd(8.116)(8.117)

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

Tracer advection schemes

Therefore,∆µ2L 
1
1
T˜ln = c00 − ξln ∆t c10 +  (ξln ∆t) 2 −
c20
2
12 
 3

1
1
1
− η ln ∆t c01 + (η ln ∆t) 2 c02 + ξln η ln ∆t 2 c11,233
∆ψ 2L 1
1
T˜dn = d 00 − η dn ∆t d 01 +  (η dn ∆t) 2 −
d 02
2
12 
 3


(8.118)1
1
1− ξ dn ∆t d 10 + (ξ dn ∆t) 2 d 20 + ξ dn η dn ∆t 2 d 11 .
2
3
3

(8.119)Finally, we describe how to derive the boundary conditions. Since the face values of the tracers are calculated throughthe second order interpolation, the value of a tracer at a point over land is sometimes necessary. For that case, the value
should be appropriately decided by using the tracer values at the neighboring points in the sea. Since ocean models
generally assume that there is no flux of tracers across land-sea boundary, the provisional value over land should be given
so as not to create a normal gradient at the boundary.
When the face value of a tracer at boundary l is calculated, W and L are not land, but either N or S may be land, andeither NW or SW may be land. When N or S is land, the land-sea boundary runs at the center of L in the zonal direction.
It is reasonable to assume that the value of land grid N or S must not cause any meridional tracer gradient at L set bysecond order interpolation using the values at grids N, L, and S. Thus, we set(TNn − TLn )∆ψ d2 = (TSn − TLn )∆ψu2 .

(8.120)

When NW or SW is a land grid, the following should be assumed.
n
n
n
n
(TNW
− TW
)∆ψ d2 = (TSW
− TW
)∆ψu2

(8.121)

n
n
n
(TWW
− TW
)∆µ2l = (TLn − TW)∆µ2ll .

(8.122)n(TEn − TLn )∆µ2l = (TW− TLn )∆µ2r .

(8.123)

When WW is a land grid,
When E is a land grid,
Similar boundary conditions are specified for face d.

8.6 Second Order Moment (SOM) scheme
8.6.1 Outline
The Second Order Moment (SOM) advection scheme by Prather (1986) seeks to improve the accuracy by treating thetracer distribution within a grid cell, unlike the scheme that aims to calculate the tracer flux at the boundary of grid cellswith high accuracy. It is assumed that the distribution of tracer f in a grid cell (0 ≤ x ≤ X, 0 ≤ y ≤ Y, 0 ≤ z ≤ Z;
volume V = XY Z) can be represented using second order functions as follows:f (x, y, z) = a0 + a x x + a xx x 2 + ay y + ayy y 2 + az z + azz z 2 + a xy xy + ayz yz + azx zx.

(8.124)

Prather (1986) expressed the above as a sum of orthogonal functions Ki (x, y, z);f (x, y, z) = m0 K0 + m x K x + m xx K xx + my Ky + myy Kyy + mz Kz + mzz Kzz + m xy K xy + myz Kyz + mzx Kzx,

– 97 –(8.125)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.6 Second Order Moment (SOM) scheme
where the orthogonal functions are given as follows:
K0 = 1,K x (x) = x − X/2,
K xx (x) = x 2 − X x + X 2 /6,
Ky (y) = y − Y /2,
Kyy (y) = y 2 − Y y + Y 2 /6,
Kz (z) = z − Z/2,

(8.126)Kzz (z) = z − Z z + Z /6,
K xy (x, y) = (x − X/2)(y − Y /2),
Kyz (y, z) = (y − Y /2)(z − Z/2),
Kzx (z, x) = (z − Z/2)(x − X/2),
2

2

Ki K j dV = 0(i , j).

∫andThe constants for normalization are decided using∫
K x2 dV = V X 2 /12,∫

(8.127)∫2
K xx
dV = V X 4 /180,∫
Ky2 dV = VY 2 /12,2
Kyy
dV = VY 4 /180,∫∫
Kz2 dV = V Z 2 /12,∫2Kzz
dV = V Z 4 /180,

∫
2
K xy
dV = V X 2Y 2 /144,∫
2
Kyz
dV = VY 2 Z 2 /144,

2Kzx
dV = V Z 2 X 2 /144.

The moments are set by the following expressions:∫S0 =f (x, y, z)K0 dV = m0V,∫
Sx = (6/X )f (x, y, z)K x (x)dV = m x V X/2,∫
Sxx = (30/X 2 )f (x, y, z)K xx (x)dV = m xx V X 2 /6,∫
Sy = (6/Y )f (x, y, z)Ky (y)dV = my VY /2,∫
Syy = (30/Y 2 )f (x, y, z)Kyy (y)dV = myy VY 2 /6,∫
Sz = (6/Z )f (x, y, z)Kz (z)dV = mz V Z/2,∫
Szz = (30/Z 2 )f (x, y, z)Kzz (x)dV = mzz V Z 2 /6,∫
Sxy = (36/XY )f (x, y, z)K xy (x, y)dV = m xy V XY /4,∫
Syz = (36/Y Z )f (x, y, z)Kyz (y, z)dV = myz VY Z/4,∫
Szx = (36/Z X )f (x, y, z)Kzx (z, x)dV = mzx V Z X/4.

(8.128)All these moments are transported with the upstream advection scheme. The procedure is carried out in one direction
at a time. The second and third procedures use the results of the last procedure. For simplicity, we describe the change– 98 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 8

Tracer advection schemesof each moment caused by an advection procedure in one direction (x) in a two dimensional plane (x y) in the following.
You may replace (y, Y ) with (z, Z ).
When velocity c in the x direction is positive, the right part of the grid cell,X − ct ≤ x ≤ X,

0 ≤ y ≤ Y,

0 ≤ z ≤ Z,

(8.129)is removed from the cell and added to the adjacent cell on the right during time interval t. This part is expressed using
superscript R. The remaining part,0 ≤ x ≤ X − ct,

0 ≤ y ≤ Y,

0 ≤ z ≤ Z,

(8.130)is expressed by superscript L. New orthogonal functions KiR (KiL ) are calculated in the part R (L) with the volume
V R = ctY Z (V L = (X − ct)Y Z). The orthogonal functions are given as follows:
K0L = K0R = 1,
K xL = x − (X − ct)/2,

K xR = x − (2X − ct)/2,

L
K xx= x 2 − (X − ct)x + (X − ct) 2 /6,
R
K xx= x 2 − (2X − ct)x + (X − ct)X + (ct) 2 /6,

KyL = KyR = y − Y /2,(8.131)

L
R
Kyy
= Kyy
= y 2 − Y y + Y 2 /6,
LK xy= [x − (X − ct)/2](y − Y /2),
RK xy
= [x − (2X − ct)/2](y − Y /2).

The basic quantities for calculating the moments areR
m0R = m0 + K̄ xR m x + K̄ xx
m xx,

m xR = m x + 2K̄ xR m xx,
R
m xx= m xx,(8.132)

myR = my + 2K̄ xR m xy,
Rmyy= myy,
Rm xy= m xy,where K̄ is the average of the new orthogonal function:
K̄ xL = −ct/2,
LK̄ xx= ct(2ct − X )/6,

K̄ xR = (X − ct)/2,
RK̄ xx= (X − ct)(X − 2ct)/6.

The moments in the right part to be removed are expressed as follows:
S0R = α[S0 + (1 − α)Sx + (1 − α)(1 − 2α)Sxx ],SxR = α 2 [Sx + 3(1 − α)Sxx ],RSxx= α 3 Sxx,

(8.133)

SyR = α[Sy + (1 − α)Sxy ],R
Syy= αSyy,R
Sxy= α 2 Sxy,– 99 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.6 Second Order Moment (SOM) scheme
where α = α R = ct/X = V R /V . The moments in the remaining part are expressed as follows:
S0L = (1 − α)[S0 − αSx − α(1 − 2α)Sxx ],SxL = (1 − α) 2 (Sx − 3αSxx ),LSxx= (1 − α) 3 Sxx,

(8.134)

SyL = (1 − α)(Sy − αSxy ),LSyy= (1 − α)Syy,
L
Sxy= (1 − α) 2 Sxy .

As the final step of the procedure, the orthogonal functions and moments transported from the adjacent cell and those in
the remaining part of the original cell are combined to create new united moments in the cell. The calculation is terribly
complex, and only the results are presented:
S0 = S0R + S0L,
Sx = αSxR + (1 − α)SxL + 3[(1 − α)S0R − αS0L ],L
R+ (1 − α) 2 Sxx+ 5{α(1 − α)(SxR − SxL ) + (1 − 2α)[(1 − α)S0R − αS0L ]},Sxx = α 2 Sxx

Sy = SyR + SyL,
R
L
Syy = Syy
+ Syy
,
R
L
Sxy = αSxy
+ (1 − α)Sxy+ 3[(1 − α)SyR − αSyL ],whereα = α R = V R /(V R + V L ).

(8.135)

The globally integrated tracer is conserved through these operations.

8.6.2 Flux limiter
Some limiters are necessary to guarantee that a tracer is positive (negative) definite.
a. Prather (1986)

Prather (1986) proposed to set limits for the moments related to the direction of advection. For instance, when the moments
are advected in the x direction,
S0 ≥ 0,Sx′ = min[+1.5S0, max(−1.5S0, Sx )],(8.136)′
Sxx= min[2S0 − |Sx′ |/3, max(|Sx′ | − S0, Sxx )],′Sxy= min[+S0, max(−S0, Sxy )].

It should be noted that the application of this limiter does not completely guarantee that the tracer will be positive (negative)
definite.
b. Merryfield and Holloway (2003)

Flux limiters introduced by Morales Maqueda and Holloway (2006) (hereinafter MH06) extended that of Prather (1986)(hereinafter PR86).
Let us consider the advection in the x direction. As in PR86, MH06 consider the mean tracer distribution (τ̄(x)) in thex direction by integrating in the y and z direction within the grid cell (Equation (2.15) of MH06).
∫ Y∫ Z1dy dz τ(x, y, z)
YZ 00[( x )2]
1x=
S0 − Sx + Sxx + 2(Sx − 3Sxx ) + 6Sxx,
VXX

τ̄(x) =

– 100 –

(8.137)

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

Tracer advection schemes

where x, y, z represent the position relative to the lower south-western corner of the grid cell and X, Y , and Z are grid
widths in the x, y, and z directions, respectively. V is the volume of the grid cell.
First we consider the minimum value of τ̄(x). When Sxx is negative (region I), the minimum is taken either at x = 0 or
x = X because τ̄(x) is convex upward. Even when Sxx is positive, if the value of x that gives the global minimum of thequadratic function is not within 0 ≤ x ≤ X (region II), the minimum is taken either at x = 0 or x = X. Otherwise (regionIII), the minimum value is given as the global minimum of the quadratic function. These are summarized as follows:
S0 − |Sx | + Sxx,
 S − |S | + S ,0xxx
V min τ̄(x) =  S − S x2 +3S x2 x ,6S x x
 0

if Sxx ≤ 0(region I),
if |Sx | ≥ 3Sxx ≥ 0 (region II),if 3Sxx ≥ |Sx |(8.138)(region III).

The line of V min τ̄(x) = 0 is plotted in the (|Sx |, Sxx ) space in Figure 8.5. If a pair of (|Sx |, Sxx ) resides in the region
bounded by the blue line and the left vertical axis, the minimum value of τ̄(x) is positive.
2S0

S0 - (S2x + 3S2xx)/(6Sxx) = 0

Sxx

S0

0.5S0

0

S0 - |Sx| + Sxx = 0
-S0
0

S0

1.5S0 1.732S0|Sx|

Figure8.5 Schematic to explain the flux limiter to avoid undershooting. This is a reproduction of Figure 3 ofMorales Maqueda and Holloway (2006). See text for explanation.

PR86 proposed to modify the moments so that the minimum of τ̄ is positive. This is to improve the sign-definiteness ofthe tracer. The limiters of MH06 are also based on this strategy. They proposed to modify the moments in the following
way,
[
[
]]
Sx′ = min 31/2 S0, max −31/2 S0, Sx ,
[
]
(
) 1/2
[ ′
]
S′2

2
x

min S0 + S0 − 3, max Sx − S0, Sxx ,if Sx′ < 23 S0,



′

[
[
]
]
(8.139)Sxx = ()() 1/2
′2 1/2S′2 min S0 + S02 − S3x, max S0 − S02 − 3x
, Sxx , if Sx′ ≥ 23 S0,[]′
Sxy = min S0, max[−S0, Sxy ] ,
[
]
′
Sxz
= min S0, max[−S0, Sxz ] .
′ , S ′ are modified moments and it is assumed that S ≥ 0. This is slightly diﬀerent from PR86 because an
Here, Sx′ , Sxx0xy√
approximation 3 ∼ 1.5 (see 8.136) was used in PR86.
This modification is explained using Figure 8.5. If a pair of (|Sx |, Sxx ) resides in the region that gives negative minimumof τ̄(x), |Sx | is shifted along the horizontal axis so that it is less than 31/2 S0 . Then Sxx is shifted along the vertical axis sothat the point is within the region that gives the minimum of τ̄(x) to be positive. The modification to Sxy, Sxz is designed
so that a negative value of the tracer is avoided by these moments.

– 101 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.6 Second Order Moment (SOM) scheme
MH06 proposed to extend the above modifications for setting the lower and upper bounds for the tracers to improve themonotonicity. The lower and upper limit for τ(x) are set as τa and τb , respectively (τa ≤ τ̄(x) ≤ τb ). The setting of thelower limit is achieved by replacing S0 with S0∗ = S0 − V τa in the above modification. For the upper limit, we interpretthe condition "the maximum of τ̄(x) must be less than τb " as "the minimum of −τ̄(x) must be more than −τb ." In otherwords, the condition is "the minimum of τb − τ̄(x) must be positive."
Using (8.137), τb − τ̄(x) may be expressed as follows:[( x )2]1xτb − τ̄(x) =V τb − S0 + Sx − Sxx − 2(Sx − 3Sxx ) − 6Sxx
VXX[( x )2]x
1 ∗S + Sx − Sxx − 2(Sx − 3Sxx ) − 6Sxx,
(8.140)=
V 0
X
X
where S0∗ = V τb − S0 . Its minimum is obtained as follows:S0∗ − |Sx | − Sxx,
S0∗ − |Sx | − Sxx,
V min [τb − τ̄(x)] =  S ∗ + S x2 +3S x2 x ,6S x x
 0

if Sxx ≥ 0(region I),
if |Sx | ≥ −3Sxx ≥ 0 (region II),if − 3Sxx ≥ |Sx |(8.141)(region III).

Following Figure 8.5, this is visualized as Figure 8.6.
S0

S*0 - |Sx| - Sxx = 0

Sxx0-0.5S0-S0S*0 + (S2x + 3S2xx)/(6Sxx) = 0

-2S0

S0

0

1.5S0 1.732S0|Sx|

Figure8.6 Schematic to explain the flux limiter to avoid overshooting. The region in (|Sx |, Sxx ) space that may set
the upper bound on tracer distribution. The region bounded by the blue curve and the left vertical axis gives the tracer
value that does not exceed the upper bound.

The modification for the moments corresponding to (8.139) is expressed as follows:[[]]
Sx′ = min 31/2 S0∗, max −31/2 S0∗, Sx ,
[
[
]]
(
)
′2 1/2
∗ − S ′ , max −S ∗ − S ∗2 − S x

min
S
,
S
,

xx
x

3
0
0
0

′

[
[
]]
Sxx
=
(
)
(
) 1/2
′2 1/2S′2 min −S0∗ + S0∗2 − S3x, max −S0∗ − S0∗2 − 3x
, Sxx ,
[]′
Sxy = min S0∗, max(−S0∗, Sxy ) ,
[
]
′
Sxz
= min S0∗, max(−S0∗, Sxz ) .

if Sx′ < 23 S0∗,
if Sx′ ≥ 23 S0∗,

(8.142)

MH06 proposed three modification methods (Method A - C). Following COCO, MRI.COM adopts "Method B" of
Morales Maqueda and Holloway (2006) as proposed by Merryfield and Holloway (2003), which lays emphasis on– 102 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 8

Tracer advection schemesmonotonicity and suppresses numerical diﬀusion as well. In this method, the minimum and the maximum value among
the three grid cells (i − 1, i, i + 1) are set as the lower and the upper limit, respectively. Moments are adjusted using
Eqs. (8.139) and (8.142). This method does not guarantee exact monotonicity, but improves monotonicity of the zeroth
moment.

8.6.3 Calculating SOM advection in MRI.COM
It should be noted that the coordinate system is not Cartesian in ocean models. Since the coordinate system covers aspherical surface, the x direction in a grid cell is not identical to that in the adjacent cell, for instance. Thus, the exact
conservation of moments cannot be realized. In addition, a tracer-cell including solid earth (sea floor or lateral boundary)is not a cuboid, so the orthogonal functions cannot be defined precisely for such a cell. Nevertheless, the procedures
described in the previous subsection can be carried out using the volume of seawater in the non-cuboid grid cell, and thezeroth moment S0 (total amount of the tracer) is conserved.
As indicated in the expressions in the previous subsections, the volume-integrated moments (Si ) and the fractionof volume to be removed (α) are used in the SOM advection scheme. There are 10 moments for each tracer. The
fraction α is calculated using volume transports (USTARL, VSTARL, and WLWL), which are calculated in the subroutines
continuity__diagnose_horizontal and continuity__diagnose_vertical. Following Prather (1986), the procedures in three directions are executed in order, not simultaneously. By default, the procedure in the meridional direction(advec_y) is called first, the zonal direction (advec_x) next, and lastly the vertical direction (advec_z). The order ofoperations in the horizontal direction may be flipped every time step if runtime option lsomstrang is set to be true (Strang
splitting). This would improve the overall accuracy of the serially executed advection operations (Skamarock, 2006). Thechange in the tracer value caused by SOM advection is estimated in the subroutine tracer__adv and added to the variable
trcal directly.
Usage

Model option SOMADVEC must be specified for compilation when the SOM advection scheme will be used for any tracer.
Namelist nml_somadv is required (see Table 8.2) at run time. You must also specify which tracer will use this scheme aswell as what kind of specifications will be used for that tracer. See Section 13.3.2 and Table 13.3 for details.
Table8.2 Variables defined in namelist nml_somadv

namefile_som_in
file_som_out
nstep_som_restart_intervallsomstrang
lsommonitor∗ i={x, xx, y, yy, z, zz, xy, yz, zx},

explanationbase name of the restart file to be input
base name of the restart file to be output
interval in time steps with which restart
files are output
The order of calling x and y direction isreversed every time step
flag to monitor the conservation of the
moments

string or value{file_som_in}_i_n ∗{file_som_out}_i_n ∗default is 0 (at the end of this job).true. / .false.
.true. / .false.

n={01, 02, ..., numtrc}; tracer number

8.7 MPDATA scheme
8.7.1 Outline
This section explains the MPDATA scheme. In the MPDATA scheme (e.g., Smolarkiewicz and Margolin, 1998), advection
for T is first solved by using the original volume transport, U T , V T , W T , to obtain a temporary value (T (1) ) using the
upstream scheme. Using this temporary value, an anti-diﬀusive volume transport (U (1) , V (1) , W (1) ) is computed. Here,the superscript (1) means that it is the first approximation to the error to be subtracted. This set of transports is used tocompute a value of the next time step using the upstream scheme starting from the above temporary value.
Original dimensionless form of the anti-diﬀusive velocity given by Smolarkiewicz and Margolin (1998) is written as– 103 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
8.7 MPDATA scheme
follows (their equation 16):
u (1) δt
= |U |(1 − U ) A(1) − U V B (1) − U WC (1),δxv (1) δt
V (1) ≡
= |V |(1 − V )B (1) − VWC (1) − VU A(1),δyw (1) δt
= |W |(1 − W )C (1) − WU A(1) − WV B (1),W (1) ≡
δz
U (1) ≡

(8.143)(8.144)(8.145)where U, V, and W are the dimensionless velocity based on the original flow field and A(1) , B (1) , and C (1) are defined as
[

] (1)
δx ∂T
A ≡
,
2T ∂ x
] (1)
[
δy ∂T
,
B (1) ≡
2T ∂ y
] (1)
[
δz ∂T
(1)
.
C ≡
2T ∂z
(1)

(8.146)(8.147)(8.148)In MRI.COM, we formulate them by using the set of original volume transports U T , V T , W T . Considering the definitionof the volume transport (8.7) to (8.9), the dimensionless form of volume transports can be written as
Ui+ 1 , j,k− 1 =
2

2

(

2

2

∆t
),

(8.149)(),∆x i, j+ 1 ∆yi, j+ 1 ∆zi− 1 , j+ 1 ,k− 1 + ∆zi+ 1 , j+ 1 ,k− 1(8.150)

2

2

2

2

2

2

2

2V T 1 1 ∆ti, j+ 2 ,k− 2

2

Wi, j,k =i+ 2 , j,k− 21∆x i+ 1 j ∆yi+ 1 , j ∆zi+ 1 , j− 1 ,k− 1 + ∆zi+ 1 , j+ 1 ,k− 1
2

Vi, j+ 1 ,k− 1 =2U T 1

2

2

Wi,Tj,k ∆t∆zi, j,k (areat)i, j,k+ 1

2

2

2

2

2

.

(8.151)

2

The A(1) through C (1) terms at (i + 21 , j, k − 21 ), which are needed to calculate U (1)1

i+ 2 , j,k− 12

A(1)1=
i+ 2 , j,k− 21∆x i+ 1 , j δ x T (1)1

i+ 2 , j,k− 21x2

2T (1)1

i+ 2 , j,k− 21

C (1)1i+ 2 , j,k− 21=

2T (1)1
i+ 2 , j,k− 12=xy∆yi+ 1 , j δ y T (1)12

i+ 2 , jxy

2T (1)1
i+ 2 , j,k− 12(8.153)+ϵxz∆zi+ 1 , j δz T (1)12

(8.152)+ϵi+ 2 , j,k− 21

B (1)1, are

i+ 2 , jxz(8.154)+ϵ

where ϵ is a small number that prevents zero division when tracer values reach zero. The finite diﬀerence and averaging
operators are defined as follows (definitions in y and z directions are the same):δ x Ai ≡x

Ai ≡Ai+ 1 − Ai− 1
2

2

∆x i
Ai+ 1 + Ai− 1
2

2

2

Ai+1 − Ai,∆x i+ 1,

δ x Ai+ 1 ≡,Ai+1 + Ai
xAi+ 1 ≡
.
2
2

2

2

The A(1) through C (1) terms should be estimated at (i, j + 21 , k − 12 ) for V (1) 1

i, j+ 2 ,k− 12– 104 –and (i, j, k) for Wi,(1)j,k .

(8.155)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 8

Tracer advection schemes

Then anti-diﬀusive volume transport is written as
U (1)1i+ 2 , j,k− 12

V (1) 1i, j+ 2 ,k− 21T= Ui+(1 − Ui+ 1 , j,k− 1 ) A(1)11, j,k− 1
2

2

2

2

i+ 2 , j,k− 12= Vi,Tj+ 1 ,k− 1 (1 − Vi, j+ 1 ,k− 1 )B (1) 1
2

2

2

2

i, j+ 2 ,k− 12

T
− Ui+
V 1
1
1
, j,k− 1 i+ , j,k−2

2

2

Wi,(1)j,k = Wi,Tj,k (1 − Wi, j,k )Ci,(1)j,k − Wi,Tj,k Ui+ 1 , j,k

xz

2

2

2

B (1)1i+ 2 , j,k− 21

2− Vi,Tj+ 1 ,k− 1 Ui, j+ 1 ,k− 1
2

xy

xy

A(1) 1i, j+ 2 ,k− 12

2

Ai,(1)j,k − Wi,Tj,k Vi+ 1 , j,k

yz

2

T
− Ui+
Wi+ 1 , j,k− 1
1, j,k− 1
2

2

2

xz

2

C (1)1i+ 2 , j,k− 12(8.156)− Vi,Tj+ 1 ,k− 1 Wi, j+ 1 ,k− 1
2

2

2

2

yz

C (1) 1i, j+ 2 ,k− 21(8.157)Bi,(1)j,k .

(8.158)8.7.2 Gauge transformation and flux limiter
Original MPDATA scheme is positive definite and cannot be used for tracers that take negative value. This can be avoided
by adding some constant for the temporary value, T (1) .
The monotonicity of MPDATA scheme can be obtained by using a non-oscillatory option proposed by Smolarkiewicz
and Grabowski (1990). In this option, an anti-diﬀusive volume transport has the upper limit given as follows:[
U (1)1] mon[
V (1) 1

] moni+ 2 , j,k− 12

i, j+ 2 ,k− 12[min(1, β out 1 , β in
)U (1)1
i+1, j,k− 12i, j,k− 2
i+ 2 , j,k− 21

=
out
in
 min(1, β,β)U (1)1
i, j,k− 12
i+1, j,k− 12
i+ 2 , j,k− 21if)V (1)
min(1, β out 1 , β ini, j+1,k− 21 i, j+ 12 ,k− 12i, j,k− 2= min(1, β out, β in 1 )V (1)1
i, j,k− 2 i+ 2 , j,k− 12
i, j+1,k− 12U (1)1> 0,

ifi+ 2 , j,k− 12
U (1)1< 0,
i+ 2 , j,k− 12

ifV (1) 1> 0,

if

i, j+ 2 ,k− 12
V (1) 1 1 < 0,i, j+ 2 ,k− 2

 min(1, β out

ifWi,(1)j,k > 0,

ifWi,(1)j,k < 0,] mon 
Wi,(1)j,k=, β in 1 )Wi,(1)j,ki, j,k− 2i, j,k+ 21
min(1, β out 1 , β in 1 )Wi,(1)j,ki, j,k+ 2
i, j,k− 2(8.159)(8.160)(8.161)where
βi,inj,k− 1 ≡2

T max 1 − T (1)i, j,k− 2i, j,k− 21,(8.162)− T min 1
i, j,k− 2i, j,k− 12,
Aout 1 + ϵ
i, j,k− 2(8.163)

Ain 1 + ϵi, j,k− 2T (1)≡βi,outj,k− 12

T max 1 = max(i, j,k− 2

T mini, j,k− 12

= min

(

Ti, j,k− 1 , Ti+1, j,k− 1 , Ti−1, j,k− 1 , Ti, j+1,k− 1 , Ti, j−1,k− 1 , Ti, j,k+ 1 , Ti, j,k− 3 ,
2
2
2
2
2
2
2
)
T (1) 1 , T (1)
T (1)
T (1)
T (1)
T (1) 1 , T (1) 3 ,
1,
1,
1,
1,

(8.164)Ti, j,k− 1 , Ti+1, j,k− 1 , Ti−1, j,k− 1 , Ti, j+1,k− 1 , Ti, j−1,k− 1 , Ti, j,k+ 1 , Ti, j,k− 3 ,
2
2
2
2
2
2
2
)
T (1) 1 , T (1)
T (1)
T (1)
T (1)
T (1) 1 , T (1) 3 ,
1,
1,
1,
1,

(8.165)i, j,k− 2i, j,k− 2

i+1, j,k− 2

i+1, j,k− 2i−1, j,k− 2i−1, j,k− 2

i, j+1,k− 2

i, j+1,k− 2i, j−1,k− 2i, j−1,k− 2i, j,k+ 2

i, j,k+ 2

i, j,k− 2i, j,k− 2

and Ain and Aout are the absolute values of the total incoming and outgoing advection flux at the T-box.

8.7.3 Usage
Model option MPDATAADVEC must be specified for compilation when MPDATA will be used for any tracer. At run time,
you must specify which tracer will use this scheme as well as what kind of specifications will be used for that tracer. See
Section 13.3.2 and Table 13.3 for details.

– 105 –,,

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 9

SGS parameterization of lateral mixing of tracersThis chapter explains subgrid-scale parameterizations for horizontal mixing of tracers.

9.1 Introduction and Formulation
Historically, a harmonic diﬀusion operator is applied in each direction of the model coordinates to express mixing oftracers. In the real ocean, transport and mixing would occur dominantly along neutral (isopycnal) surfaces. Thus,
horizontal mixing along a constant depth surface is generally inappropriate since neutral surfaces are generally slanting
relative to a constant depth surface. Neutral physics schemes are devised as substitutes for the harmonic scheme in thehorizontal direction, while the harmonic scheme continues to be used for vertical diﬀusion.
By default, the diﬀusion operator mixes a tracer in each direction of the model coordinates with the harmonic scheme.
For horizontal diﬀusion, the biharmonic scheme can be used instead of the harmonic scheme. Using (20.22) and (20.23),
the harmonic-type diﬀusivity is represented for any tracer (T) as follows:
{
(
)
(
)}
1
∂ hψ z s κ H ∂T
∂ hµ z s κ H ∂T
∂ ( κV ∂T )−z s ∇ · FT = D(T ) =++h µ hψ ∂ µhµ ∂ µ∂ψhψ ∂ψ
∂s z s ∂sTTT ∂(hψ z s Fµ ) ∂(hµ z s Fψ ) 
 ∂Fz
1 

−
=−
+
,
(9.1)


h µ hψ
∂µ∂ψ∂swhere κ H and κV are the horizontal and vertical diﬀusion coeﬃcients and the diﬀusive fluxes are represented by F. D(T )represents the diﬀusion operator for T. When the biharmonic-type is selected for horizontal diﬀusion, the above Laplacian
operation is repeated twice reversing its sign after the first operation.
When the neutral physics schemes are selected, the advection-diﬀusion equation for any tracer (T) is expressed as follows(Gent and McWilliams, 1990):
[ ∂
]
]
DT
∂ [
+ ∇ H · T (κT ∇ H ρ/ρz ) +T∇ H · (−κT ∇ H ρ/ρz ) = D(T ) + QT ,
Dt
∂z
∂z

(9.2)

where the first term on the r.h.s. is the isopycnal diﬀusion, whose form is given byD(T ) = ∇ · (κ I K∇T ),

(9.3)

where
ρ2 + ρ2z
*. y
1
K= 2
. −ρ x ρyρ x + ρ2y + ρ2z −ρ ρx z,

−ρ x ρyρ2x + ρ2z
−ρy ρz−ρ x ρz+−ρy ρz //ρ2x + ρ2y -1 + ( ρy /ρz ) 2
1*. −( ρ x /ρz )( ρy /ρz )=1 + ( ρ x /ρz ) 2 + ( ρy /ρz ) 2−ρ x /ρz,(9.4)−( ρ x /ρz )( ρy /ρz )1 + ( ρ x /ρz ) 2−ρy /ρz−ρ x /ρz+/ ,−ρy /ρz22( ρ x /ρz ) + ( ρy /ρz ) -(9.5)(Redi, 1982). In the above, the Cartesian notation is used for simplicity. The subscript x represents ∂/(hµ ∂ µ) y represents∂/(hψ ∂ψ), and z represents ∂/(z s ∂s), The isopycnal diﬀusion coeﬃcient is represented by κ I . Diapycnal diﬀusion is notconsidered here.

– 107 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
9.2 Horizontal diﬀusion
The second and third terms on the l.h.s. of (9.2) have the form of advection terms with a transport velocity vector(uT , vT , wT ):∂ (1 ∂ρ ∂ρ)/κT,∂zhµ ∂ µ ∂z
1 ∂ρ ∂ρ)
∂ (
/
κT
,
vT ≡
∂z
hψ ∂ψ ∂z
1 { ∂ ( hψ ∂ ρ ∂ ρ )∂ ( hµ ∂ ρ ∂ ρ )}wT ≡ −
/
/
κT
+
κT
,
h µ hψ ∂ µhµ ∂ µ ∂z∂ψhψ ∂ψ ∂z
uT ≡

(9.6)
(9.7)
(9.8)

(Gent and McWilliams, 1990). This velocity can be understood as the advection caused by the thickness diﬀusion of an
isopycnal layer, where κT is often referred to as thickness diﬀusivity.
Note that these could be rewritten asG(T ) = ∇ · (κT A∇T )
(9.9)
with
0
A = *. 0, ρ x /ρz00ρy /ρz−ρ x /ρz−ρy /ρz +/ .
0-

(9.10)

Comparing with (9.5), we notice that the isopycnal diﬀusion and the thickness diﬀusion terms are combined to yield aconcise form (Griﬃes, 1998) and (9.2) can be rewritten as:DT= ∇ · {(κ I K − κT A)∇T } + QT .
Dt

(9.11)

Three types of horizontal diﬀusion, harmonic horizontal diﬀusion (default), biharmonic horizontal diﬀusion (TRCBIHARM
option), and isopycnal diﬀusion (ISOPYCNAL option), are available in MRI.COM. When isopycnal diﬀusion (Redi, 1982)
is selected, the parameterization scheme for eddy induced advection by Gent and McWilliams (1990) (GM scheme) isused with it (see Section 9.4 for details of this parameterization).
The following is a guide to selecting a horizontal diﬀusion scheme. Biharmonic diﬀusion is appropriate for a high resolution model that can resolve eddies because it is more scale-selective than harmonic diﬀusion and does not unnecessarily
suppress disturbances in resolved scales. On the other hand, biharmonic diﬀusion is not recommended in eddy-less models
because this would result in numerical instability. Harmonic horizontal diﬀusion is also not recommended because thisscheme would cause unrealistic cross-isopycnal (diapycnal) mixing as mentioned above. Instead, both isopycnal diﬀusion
and the GM scheme should be used there. Using an anisotropic GM scheme can maintain the meso-scale eddy structuresand swift currents by restricting the direction of diﬀusion, and thus may be usable even for a high resolution model.
The finite diﬀerence expression of (9.11) is given by taking the finite volume approach as follows:
Ti,n+1∆Vi,n+1= Ti,n−1∆Vi,n−1+ 2 ∆t {FXDi− 1 , j,k− 1 − FXDi+ 1 , j,k− 1 + FYDi, j− 1 ,k− 1 − FYDi, j+ 1 ,k− 1j,k− 1j,k− 1j,k− 1j,k− 1
2

2

2

2

2

2

2

2

2

+ FZDi, j,k − FZDi, j,k−1 } + (other terms).

2

2

2

(9.12)

Usually, fluxes due to diﬀusion are explicitly represented using a starting time level of the temporal discretization. However,when the flux is very large relative to the grid size and the time step interval, which would often occur for vertical fluxes,
an implicit scheme is used. These issues are discussed in Chapter 19.

9.2 Horizontal diﬀusion
9.2.1 Laplacian diﬀusion
Harmonic horizontal diﬀusion assumes that the diﬀusion flux is a product of the gradient of tracer and the diﬀusion
coeﬃcient (κ H ). The finite diﬀerence expressions for the fluxes are given as follows:y

FXDi+ 1 , j,k− 1 = − κ H ∆yi+ 1 , j ∆zi+ 1 , j,k− 1 δ x Ti+ 1 , j,k− 1 ,
2

2

2

2

2

2

2

x

FYDi, j+ 1 ,k− 1 = − κ H ∆x i, j+ 1 ∆zi, j+ 1 ,k− 1 δ y Ti, j+ 1 ,k− 1 ,2

2

2– 108 –

2

2

2

2

(9.13)(9.14)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 9

SGS parameterization of lateral mixing of tracers

where
Ti+1, j,k− 1 − Ti, j,k− 1

δ x Ti+ 1 , j,k− 1 ≡22

2∆x i+ 1 , j

2,(9.15).

(9.16)

2

Ti, j+1,k− 1 − Ti, j,k− 1

δ y Ti, j+ 1 ,k− 1 ≡222∆yi, j+ 1

2

2

9.2.2 Biharmonic diﬀusion
Biharmonic horizontal diﬀusion (TRCBIHARM option) assumes that the diﬀusion flux is proportional to the gradient of the
Laplacian of tracer. The finite diﬀerence expressions for the fluxes are given as follows:√y(9.17)
FXDi+ 1 , j,k− 1 = κ b ∆yi+ 1 , j ∆zi+ 1 , j,k− 1 δ x [Lap (T )]i+ 1 , j,k− 1 ,
2
2
2
2
2
2
2
√
x
FYDi, j+ 1 ,k− 1 = κ b ∆x i, j+ 1 ∆zi, j+ 1 ,k− 1 δ y [Lap (T )]i, j+ 1 ,k− 1 ,(9.18)
2

2

2

2

2

2

2

where κ b is diﬀusivity coeﬃcient andLap (T )i, j,k− 1 =2∆x i, j ∆yi, j√√yx(δ x κ b ∆zi, j,k− 1 δ x Ti, j,k− 1 + δ y κ b ∆zi, j,k− 1 δ y Ti, j,k− 1 ).
2
2
2
2
∆Vi, j,k− 1(9.19)

2

9.2.3 Specification of coeﬃcient
The diﬀusion coeﬃcient of horizontal diﬀusion is specified using the namelists listed on Tables 9.1 and 9.2. This must be
zero if ISOPYCNAL option is selected, unless the horizontal diﬀusion is applied intentionally.
Table9.1 Namelist nml_tracer_diff_horz

variable name
diff_horz_cm2ps
diff_horz_cm4ps
file_diff_horz_2d

units
cm2 s−1
cm4 s−1characterdescription
Laplacian diﬀusion coeﬃcient (κ H )Biharmonic diﬀusion coeﬃcient (κ b )2D distribution of diﬀusion coeﬃcient is read fromthis fileusage
if not TRCBIHARM
if TRCBIHARM

Table9.2 Namelist nml_grid_size_change_mix_coefs

variable namel_grid_size_change_mix_coefs

units
logicaldescription
the given coeﬃcient is multiplied by thefraction of the local grid size to 100 km.

usage

9.3 Isopycnal diﬀusion
In isopycnal diﬀusion, the diﬀusion flux is expressed by high diﬀusivity along the isopycnal surface κ I , low diﬀusivity
across the isopycnal surface κ D , and their product with the gradient of tracer in each direction. Using diﬀusion tensor K,
each flux component is written as F m (T ) = −K mn ∂nT, and then
1 + Sy2 + ϵ Sx2*.
K=. (ϵ − 1)Sx Sy
1 + S2, (1 − ϵ )SxκI(ϵ − 1)Sx Sy
1 + Sx2 + ϵ Sy2(1 − ϵ )Sy(1 − ϵ )Sx+(1 − ϵ )Sy // ,ϵ + S2 -(9.20)

where
∂ρ ∂ρ
*.
+
∂ y //, 0/ ,
S = (Sx, Sy, 0) = ..− ∂ x , −. ∂ρ ∂ρ /S = | S |,
κD
ϵ=
κI

, ∂z

∂z

(9.21)(9.22)(9.23)– 109 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
9.3 Isopycnal diﬀusion(Redi, 1982).
The finite diﬀerence method is based on Cox (1987) except for the small isopycnal slope approximation. The finite
diﬀerence form of three components of the gradient of tracer in calculating the east-west component of flux FXDi+ 1 , j,k− 1
2
2
is given as follows (see Figure 9.1):(δ x T )i+ 1 , j,k− 1 = δ x Ti+ 1 , j,k− 1 ,2

2

2(9.24)2(δ y T )i+ 1 , j,k− 1 = δ y Ti+ 1 , j,k− 1xy,(9.25)(δz T )i+ 1 , j,k− 1 = δz Ti+ 1 , j,k− 1xz.

(9.26)xy,(9.27)

2

2

2

2

2

2

2

2

Similarly, the north-south component FYDi, j+ 1 ,k− 1 is given by22(δ x T )i, j+ 1 ,k− 1 = δ x Ti, j+ 1 ,k− 1
2

2

2

2

(δ y T )i, j+ 1 ,k− 1 = δ y Ti, j+ 1 ,k− 1 ,
2

2

2

2

yz

(δz T )i, j+ 1 ,k− 1 = δz Ti, j+ 1 ,k− 1 ,
2

2

2

2

(9.28)(9.29)and the vertical component FZDi, j,k is given byxz(δ x T )i, j,k = δ x Ti, j,k ,(9.30)yz(δ y T )i, j,k = δ y Ti, j,k ,(δz T )i, j,k = δz Ti, j,k .

(9.31)(9.32)Ti, j+1

Ti+1, j+1

Ti, j

Ti+1, jTi, j-1

Ti+1, j-1

Figure9.1 The way of calculating the gradient at the circle (i + 12 , j, k − 21 ) in isopycnal diﬀusion: the east-west
gradient is indicated by an arrow through the circle, and the north-south gradient is given by averaging four arrows in
the vertical direction.

The density gradient in the calculation of each component of the diﬀusion tensor can be obtained by replacing T in
the above equations with ρ. However, density is calculated at the reference level k − 12 for the east-west and north-south
components, and at the reference level k for the vertical component.
In addition, the upper bound on the isopycnal slope Smax (slope_clip_iso in namelist nml_isopy_slope_clip,
Table 9.4) is set because a nearly vertical isopycnal slope and the resultant low horizontal diﬀusivity could cause numericalinstability. If |S| > Smax , ∂z ρ in all flux components is replaced so as to satisfy |S| = Smax .
– 110 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 9

SGS parameterization of lateral mixing of tracers

The vertical flux due to the third diagonal component of the diﬀusion tensor (9.20) is
FZDi, j,k = (areat)i, j,kThus the eﬀective vertical diﬀusivity κ eﬀ isκ eﬀ =

κ I (ϵ + S 2 )δz Ti, j,k .
1 + S2

κ D + κ I S2.
1 + S2(9.33)(9.34)For a steep isopycnal slope S ∼ 1/100 and a canonical value of isopycnal diﬀusion coeﬃcient κ I ∼ 107 cm2 s−1 and atypical value of diapycnal diﬀusion coeﬃcient κ D ∼ 10−1 cm2 s−1 ,κ eﬀ ∼ 103 cm2 s−1 .

(9.35)

This is a fairly large value which warrants use of an implicit scheme (Section 19.5). In MRI.COM, this term is separated
from other terms and solved with other vertical diﬀusion terms using an implicit method.
Griﬃes et al. (1998) noted a problem in the finite diﬀerence expression of the isopycnal diﬀusion as implemented in
the GFDL-model by Cox (1987). The problem is that the down-gradient orientation of the diﬀusive fluxes along theneutral directions does not necessarily guarantee the zero isoneutral diﬀusive flux of locally referenced density (e.g.,potential temperature when it is the only active tracer). This is caused by the nature of the finite diﬀerence method andthe non-linearity of the equation of state. Griﬃes et al. (1998) proposed a remedy, but this remains to be implemented in
MRI.COM.

9.3.1 Tapering around steep isopycnal slopes
As noted above, we set an upper bound on the isopycnal slope used to evaluate isopycnal tracer diﬀusion terms in MRI.COM
in order to prevent numerical instability around steep isopycnal slopes. Griﬃes (2004) shows that such slope clipping
could lead to an unrealistically large tracer flux. Danabasoglu and McWilliams (1995) propose another method that uses
a hyperbolic tangent to exponentially taper isopycnal diﬀusion in steep slope regions. The following factor is multiplied
to the isopycnal diﬀusion coeﬃcient there:{()}
Scenter − |S|1
f steep =1 + tanh.
(9.36)
2
SwidthTwo parameters Scenter and Swidth , which determine a transitional region, are given by center_transition andwidth_transition in namelist nml_tracer_diff_isopy_taper, respectively (Table 9.5). This tapering method
is valid when ISOTAPER option is chosen and used with the other factor f surface introduced in the subsequent section.

9.3.2 Surface tapering
The vertical displacement of a water parcel due to mesoscale eddy stirring (D) is approximately calculated as follows:
D = R|S|,(9.37)

where R is the internal deformation radius and |S| is the isopycnal slope. If the depth of the water parcel (d) is shallower
than D, the boundary constrains its displacement and eddy diﬀusive fluxes (Griﬃes, 2004). The ISOTAPER option enables
the use of the following tapering scheme proposed by Large et al. (1997):
{
(
)}
1
1
1 + sin π r −,(9.38)
f surface =22where r = max(0, min(1, d/D)). Eddy diﬀusivity is tapered oﬀ to zero toward the sea surface in the region where
0 ≤ d ≤ D. In MRI.COM, the reference level of the depth d is set at the boundary layer depth (BLD): d = −z − BLD. In
MRI.COM, the surface mixed layer depth (MLD) is treated as the BLD (see also the next section). This means that theeddy diﬀusivity is tapered to zero within the mixed layer. The upper boundary of this tapering region can be changed to a
constant level by setting upper_level_isotaper_m in namelist nml_tracer_diff_isopy_taper (Table 9.5).
By default, the horizontal diagonal terms of the isopycnal diﬀusion tensor (9.20) are not multiplied by the two factors
noted above. This means that the isopycnal diﬀusion is rendered the horizontal diﬀusion within the boundary layer
and in the steep slope region. One may apply these tapering methods to the horizontal diagonal terms by setting
l_apply_hdiag = .true. in namelist nml_tracer_diff_isopy_taper (Table 9.5).
– 111 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
9.4 Gent and McWilliams parameterization for eddy-induced transport

9.3.3 Specification of coeﬃcient
The diﬀusion coeﬃcients of isopycnal diﬀusion and GM parameterization explained in the next section are specified using
the namelist listed on Table 9.3. We can use diﬀerent slope maximal limits for isopycnal diﬀusion and GM parameterizationthat are specified using the namelist listed on Table 9.4. Configurations of surface tapering for the isopycnal diﬀusion
scheme with ISOTAPER option are specified using the namelist listed on Table 9.5.
Table9.3 Namelist nml_tracer_diff_isopy

variable name
diff_isopy_cm2ps
diff_diapy_cm2ps
diff_thick_cm2psunits
cm2 sec−1
cm2 sec−1
cm2 sec−1

description
isopycnal diﬀusion coeﬃcient
diapycnal diﬀusion coeﬃcient
coeﬃcient of GM parameterization

usage
if ISOPYCNAL
if ISOPYCNAL
if ISOPYCNAL

Table9.4 Namelist nml_isopy_slope_clip

variable nameslope_clip_iso

units
1

slope_clip_gm

1

description
maximum slope of isopycnal surface for isopycnal diﬀusion
maximum slope of isopycnal surface for GM
parameterization

usage
if ISOPYCNAL
if ISOPYCNAL

Table9.5 Namelist nml_tracer_diff_isopy_taper

variable namel_apply_hdiag

units
logicalcenter_transition1

width_transition1

ai_minupper_level_isotaper_m
l_explicit_vdif

cm2 sec−1mlogicaldescription
apply tapering factors to horizontal diagonal terms (default = .false.)
center of the transition region of taperingwith hyperbolic tangent (default = 0.005)
width of the transition region of taperingwith hyperbolic tangent (default = 0.001)
lower limit of horizontal isopycnal diﬀusion coeﬃcient when l_apply_hdiag =.true. (default = 0.0)
BLD for the sine taper (default = mixed
layer depth (−1.0))
handle the vertical diﬀusion term explicitly(default = .false.)

usage
if ISOTAPER
if ISOTAPER
if ISOTAPER
if ISOTAPER

if ISOTAPER
if ISOTAPER

9.4 Gent and McWilliams parameterization for eddy-induced transport
9.4.1 General features
The Gent and McWilliams (1990) parameterization represents transports of tracers due to disturbances smaller than thegrid size, assuming that a flux proportional to the gradient of the layer thickness exists along the isopycnal surface. The
isopycnal diﬀusion stated above does not produce any flux when the isopycnal surface coincides with the isotherm andisohaline surface. This parameterization, however, produces fluxes in such a case, and acts to decrease the isopycnal slope.
Flux convergence due to diﬀusion is expressed as follows:
R(T ) = ∂m (J mn ∂nT )(9.39)

Diﬀusion tensor J mn is expressed as the sum of the symmetric component K mn = (J mn + J nm )/2 and the anti-symmetriccomponent Amn = (J mn − J nm )/2. Isopycnal diﬀusion has the form of a symmetric diﬀusion tensor. Convergence of a

– 112 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 9

SGS parameterization of lateral mixing of tracersm = −Amn ∂ T is as follows:
skew flux caused by the anti-symmetric component Fskewn

R A (T ) = ∂m ( Amn ∂nT )= ∂m ( Amn )∂nT
= ∂n (∂m AmnT ),(9.40)

where Amn ∂m ∂nT = 0 and ∂m ∂n Amn = 0 are used. If we set a virtual velocity u∗n ≡ −∂m Amn , then the flux due to theanti-symmetric component could be regarded as the advection due to this virtual velocity u∗ . In this case, the flux is
Fadv = u∗T and R A (T ) = −u∗ · ∇T since u∗ is divergence free by definition.
The Gent and McWilliams parameterization for eddy-induced transport velocity is given byu∗ = −/ ∂h∂(κ GM ∇ρ h),∂ρ∂ρwhere h is the depth of the neutral surface (ρ = const). This velocity is expressed in the depth coordinate as()
u∗ = −∂z (κ GM Sx ), −∂z (κ GM Sy ), ∇h · (κ GM S)
where

(
)
S = (Sx, Sy, 0) = −ρ x /ρz , −ρy /ρz , 0 .

(9.41)(9.42)(9.43)(Gent et al., 1995).
Griﬃes (1998) showed that the tendency of a tracer due to this parameterization might be expressed using an antisymmetric diﬀusion tensor A00−κ GM Sx00−κ GM Sy +/ ,A = *.
(9.44)
κ
S
κ
S
0
GM y
, GM xso that∂T= · · · + ∇ · (A∇T ).
∂t(9.45)

The flux due to advection can be expressed using a vector streamfunction,Ψ = κ GM ẑ × S = (−κ GM Sy, κ GM Sx, 0),
which produces u∗ in (9.42):(9.46)

Fadv = u∗T = T (∇ × Ψ).

The skew diﬀusive expression for the flux using tensor A in (9.44) is
Fskew = −A∇T = −(∇T ) × Ψ = Fadv − ∇ × (T Ψ).
Thus, the convergence of the flux expressed in tensorial form matches that of the advective expression. In other words, theGent and McWilliams parameterization is realized by only adding A to the tensor of the isopycnal diﬀusion K (Griﬃes,1998).

9.4.2 Dependency of coeﬃcient on space and time
By default, the diﬀusivity coeﬃcient for the Gent-McWilliams parameterization is constant both in time and space, whose
value, [κ GM ]ref , is given by diff_thick_cm2ps in namelist nml_tracer_diff_isopy (Table 9.3). However, it may be
dependent on local horizontal grid size by specifying a namelist (see the next paragraph). Several parameterization may
be used by choosing GMVAR option. User should specify one of l_visbeck, l_eden, l_danabasoglu to be .true. in
namelist nml_gmvar_select (Table 9.8).
a. Simple scheme

If l_grid_size_change_mix_coefs = .true. in namelist nml_grid_size_change_mix_coefs (Table 9.7), the coefficient may be dependent on the horizontal grid size according to the following formulaκ GM = [κ GM ]ref × min(∆x, ∆y)/(100 km),where ∆x and ∆y are local zonal and meridional grid sizes of a U-cell, respectively.
– 113 –

(9.47)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
9.4 Gent and McWilliams parameterization for eddy-induced transport

b. Visbeck et al. (1997)

To use the method proposed by Visbeck et al. (1997), l_visbeck = .true. in namelist nml_gmvar_select.
Visbeck et al. (1997) proposed to give the GM coeﬃcient κ GM asκ GM = αwhere α = 0.015,
M2 =

M2 2
l ,
N

g
∇H ρ ,
ρ0(9.48)

N2 = −

g ∂ρ
,
ρ0 ∂z

and l the horizontal length scale of the baroclinic zone, g acceleration of gravity, ρ0 reference density.
Specifically in our model,) 2 (∫ −H0) 2  12 (∫ −H0

g
∂
ρ
∂
ρ

M2 =
dz +
dz  ,
ρ0 D  −H1 ∂ x
∂
y

−H1

and
N2 =

g [σ0 (H1 ) − σ0 (H0 )],ρ0 D(9.49)(9.50)

where H0 = 100 m, H1 = 2000 m, D = H1 − H0 , and σ0 is the potential density. Lower limit for N is set so thatN 2 ≥ 10−9 s−1 .
Using the following formula for the phase speed of the 1st baroclinic mode gravity wave (Sueyoshi and Yasuda, 2009)
1
c1 =
π

∫ 0

(

−H B

g ∂σ0
−
ρ0 ∂z

) 12
dz,

where HB is the depth of sea floor, deformation radius λ 1 is calculated as follows:()
c14
λ 1 = min, 4 × 10m.
f

(9.51)(9.52)

Using a factor r = 7, GM coeﬃcient is determined as follows:κ GM = αM2(r λ 1 ) 2 .
N(9.53)

Lower and upper limits for the coeﬃcient are set as follows:
300 ≤ κ GM ≤ 1500

m2 s−1 .

(9.54)

c. Eden and Greatbatch (2008)

To use the method proposed by Eden and Greatbatch (2008), l_eden = .true. in namelist nml_gmvar_select.
Eden and Greatbatch (2008) and Eden et al. (2009) proposed that the thickness diﬀusivity is given byκ GM = cL 2 σ(9.55)The eddy length scale L is given as the minimum of the Rossby radius L r and Rhines scale L Rhi . This choice for L wasfound to be consistent with independent estimates of eddy length scales from satellite observations and high-resolution
model results (Eden, 2007) and with recent theoretical considerations (Theiss, 2004). L Rhi is estimated from variables of
the coarse resolution model asσ(9.56)
L Rhi =
β
(Eden and Greatbatch, 2008), while L r is given by[c1,
L r = min|f|
– 114 –

√

]
c1
,
2β

(9.57)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 9

SGS parameterization of lateral mixing of tracers

where c1 denotes the 1st baroclinic gravity wave speed calculated approximately as eq. (9.51). Considering the thermal
wind relation in mid-latitudes, Eden and Greatbatch (2008) proposed that the inverse eddy time scale σ is given by
σ = f (Ri + γ) − 2 .
1(9.58)∂
Here, Ri = N 2 | ∂zuh | −2 denotes the local Richardson number. γ(> 0) is introduced to prevent the singularity as N → 0,
which acts eﬀectively as an upper limit for σ and consequently for κ GM . The default values of γ and c in eq. (9.55), are
200 and 2, respectively.

d. Danabasoglu and Marshall (2007)

To use the method proposed by Danabasoglu and Marshall (2007), l_danabasoglu = .true. in namelist
nml_gmvar_select.
Danabasoglu and Marshall (2007), guided by Ferreira et al. (2005) and Ferreira and Marshall (2006), proposed to
specify the vertical variation of κ GM using[ N2 ]
κ GM =[κ GM ]ref,(9.59)
2
Nref
where N is the local buoyancy frequency and [κ GM ]ref is the constant reference value of κ GM within the surface diabatic
layer. Nref is the reference buoyancy frequency obtained just below the diabatic layer, in other words, the first stable N 2
2 is set to 1 for all shallower depths. Between the depth at which N 2 = N 2
below surface diabatic layer. The ratio N 2 /Nref
refand the ocean bottom, we also ensure that
N2
Nmin ≤ 2 ≤ 1.0,
(9.60)
Nref
where Nmin is the lower limit specified by the user (ratio_bvf_min in namelist nml_gmvar_danabasoglu).

9.4.3 Surface tapering
By default, no specific modification is applied to the eddy-induced transport velocity of the Gent-McWilliams parameterization near the surface and the bottom, except for limiting the isopycnal slope to a specified value (slope_clip_gm innamelist nml_isopy_slope_clip). This may result in too strong transport velocity in the first vertical level of the model(sea surface). The problem may be overcome by tapering the transport in the surface mixed layer, where the transport
is made nearly or completely uniform in the vertical direction. This is realized by choosing either SLIMIT or GMTAPER
option.
a. Simple scheme

By choosing SLIMIT option, the Gent-McWilliams coeﬃcient (κ GM ) is linearly reduced from the value at the base of
mixed layer to zero at the sea surface within the mixed layer,
κ GM (z) = κ GM (z = −MLD) × (−z/MLD) for − MLD ≤ z ≤ 0,

(9.61)

where MLD is the depth of the mixed layer. The MLD is defined as the level at which the local potential density is larger
than the surface density by a specified value, given by the user (default value is 0.03 kg m−3 ).
b. Danabasoglu et al. (2008)

By choosing GMTAPER option, a practical scheme proposed by Danabasoglu et al. (2008) is used. This scheme modifies
the Gent-McWilliams vector stream function for eddy induced transport velocity near the surface, aiming to implement anear-surface parameterization proposed by Ferrari et al. (2008). Concept of the near-surface parameterization is as follows(Danabasoglu et al., 2008):• In the turbulent boundary layer (BL), the eddy-induced velocity is set parallel to the boundary and has no vertical
shear, as expected in the mixed layer.
• There is an eddy diﬀusion of buoyancy along the boundary as well as along isopycnals.
• In the interior the parameterization satisfies the adiabatic constraint as in the original scheme.
• The two forms are matched through a transition layer that separates the quasi-adiabatic interior with isopycnicallyoriented eddy fluxes from the near boundary regions.

– 115 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
9.4 Gent and McWilliams parameterization for eddy-induced transport
Two vertical length scales must be estimated to implement this parameterization: the boundary layer depth (BLD) andthe transition layer thickness (TLT). Their sum is defined as the diabatic layer depth (DLD), over which the upper-ocean
eddy fluxes depart from their interior formulas. In MRI.COM, the surface mixed layer depth (MLD) is treated as the BLD.
The MLD is defined as the level at which the local potential density is larger than the surface density by a specified value,
given by the user (default value is 0.03 kg m−3 ). The TLT is defined by the range of isopycnals that can be lifted into theboundary layer by eddy heaving, which is given by the product of the internal deformation radius (R) and the isopycnal
slope (|S|):
D = R|S|.
(9.62)
Thus we calculate D at each grid point and the DLD is obtained as follows:
DLD = BLD + D.

(9.63)

Now the near-surface expression for the eddy-induced vector streamfunction is given in the following. The streamfunctionis split into its boundary layer, ΨBL , and transition layer, ΨTL , expression as follows:
ΨBL =
and
ΨTL =

( z + BLD ) 2
TLTη−zΨ0 for − BLD ≤ z ≤ ηη + BLDΦ+( η−z )
Ψ0 for − DLD ≤ z < −BLD
η + BLD(9.64)(9.65)The two functions Ψ0 and Φ are chosen such that Ψ and its vertical derivative are continuous across the base of BLD
and the base of TLT. These constraints then yieldη + BLD(2ΨI + TLT∂z ΨI )2(η + BLD) + TLT(9.66)TLT(ΨI + (η + DLD)∂z ΨI ),2(η + BLD) + TLT(9.67)

Ψ0 =
and
Φ=−

where ΨI is the interior eddy-induced streamfunction at the base of the transition layer given by the Gent-McWilliams
parameterization,z × ∇H ρ
ΨI = −κ GM
at z = −DLD.
(9.68)∂z ρ
In the implementation, to evaluate both ΨI and ∂z ΨI at z = −DLD, ΨI are evaluated at the vertical grid points thatstraddle z = −DLD.
Danabasoglu et al. (2008) also showed that the model solutions are not very sensitive to their transition layer
thickness. Whether the transition layer is included or not may be specified by l_transition_layer in namelist
nml_gm_transition.

9.4.4 Anisotropic Gent-McWilliams scheme
An anisotropic GM scheme (Smith and Gent (2004), GMANISOTROP option), which gives greater diﬀusivity only in the
direction of the current vector, is also available. Using unit vector n̂ = (n x, ny ) in an arbitrary direction, the two-dimensional
anisotropic diﬀusion tensor is defined as follows:() ()κ A n2x + κ B n2y
κ B n x ny
L M
K2 =
=
,
(9.69)M N
κ B n x ny
κ B n2x + κ A n2ywhere κ A is the diﬀusivity in the n̂ direction, and κ B is that in the direction normal to n̂. This is applied to the antisymmetric tensor in the Gent-McWilliams scheme, and the following expression is obtained (Smith and Gent, 2004),00A′ = *.
, LSx + M Sy00
M Sx + N Sy

−LSx − M Sy
−M Sx − N Sy +/ .
0-

(9.70)In the choice of GMANISOTROP option, n̂ is set in the direction of the local horizontal velocity. The value of κ A is readfrom namelist nml_tracer_diff_isopy (variable name diff_thick_cm2ps). The ratio of κ B /κ A is read from namelist
nml_gmanisotrop (variable name cscl_isotrop). The default value of cscl_isotrop is set to 1/2.
– 116 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 9

SGS parameterization of lateral mixing of tracers

9.4.5 Usage Summary
How to specify the overall behavior of the Gent-McWilliams parameterization is summarized as follows.
a. Model options

Model options related to the GM parameterization are listed on Table 9.6
Table9.6 List of model options related to GM parameterization.

option nameGMVAR
SLIMIT
GMTAPERGMANISOTROPAFC

description
Coeﬃcient of GM parameterization is allowed to varyLinearly reduce the coeﬃcient of GM parameterizationfrom the bottom of the mixed layer to the sea surface
Taper GM vector stream function near the sea surface
An-isotropic horizontal variation of GM parameterization
Calculate additional flux by using horizontal gradients ofdensity and velocity (Hirabara et al., 2010)

usage
specify nml_gmvar
cannot be used with GMTAPER
cannot be used with SLIMIT,
GMANISOTROP, AFC
specify nml_gmanisotrop
cannot be used with TRCBIHARM

b. Spatial dependency

The diﬀusion coeﬃcient of GM parameterization may be grid size dependent by using the namelist listed on Table 9.7.
Table9.7 Namelist nml_grid_size_change_mix_coefs

variable namel_grid_size_change_mix_coefs

units
logicaldescription
the given coeﬃcient is multiplied by thefraction of the local grid size to 100 km.

usage
default = .false.

Overall behavior of GM parameterization with GMVAR option should be specified by using the namelists listed on Tables
9.8 through 9.12.
Table9.8 Namelist nml_gmvar_select for GMVAR

variable namel_visbeckl_edenl_danabasoglu

units
logicallogicallogicaldescription
use Visbeck et al. (1997)
use Eden and Greatbatch (2008)
use Danabasoglu and Marshall (2007)

usage
choose only one of the three options
choose only one of the three options
choose only one of the three options

Table9.9 Namelist nml_gmvar_visbeck for GMVAR

variable namestart_depth
base_depth
cscl_gmvar
upper_limit
lower_limit
lcalc_defrad

units
cm
cm
1
cm2 sec−1
cm2 sec−1
logicaldefrad_const

cm

length_factor

1

description
density gradients are averaged from start_depth
density gradients are averaged to base_depth
parameter for GM diﬀusivity calculation
upper limit of thickness diﬀusivity
lower limit of thickness diﬀusivity
flag whether deformation radius is calculated ornotupper limit of deformation radius whenlcalc_defrad = .true.
constant horizontal length scale whenlcalc_defrad = .false.
[horizontal length scale] = [deformation radius] ×length_factor– 117 –

usage
l_visbeck = .true.
l_visbeck = .true.
l_visbeck = .true.
l_visbeck = .true.
l_visbeck = .true.
l_visbeck = .true.
l_visbeck = .true.

l_visbeck = .true.

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
9.4 Gent and McWilliams parameterization for eddy-induced transport
Table9.10 Namelist nml_gmvar_eden for GMVAR

variable name
c_EG

units
1

gamma_EG
upper_limit
lower_limit

1
cm2 sec−1
cm2 sec−1

description
In this parameterization, Thickness diﬀusivity isparameterized as c_EG × L 2 × sigma_EG.
sigma_EG = f /(Ri + gamma_EG)
upper limit of thickness diﬀusivity
lower limit of thickness diﬀusivity

usage
l_eden = .true.
l_eden = .true.
l_eden = .true.
l_eden = .true.

Table9.11 Namelist nml_gmvar_danabasoglu for GMVAR

variable nameratio_bvf_min

units
1

ratio_bvf_max

1

description
Lower bound for the squared buoyancy frequency
relative to the reference value (default = 0.1)
Upper bound for the squared buoyancy frequency
relative to the reference value (default = 1.0)

usage
l_danabasoglu = .true.
l_danabasoglu = .true.

Table9.12 Namelist nml_gm_transition

variable namel_transition_layer

units
logicaldescription
include transition layer into diabatic
layer (default = .false.)

usage
eﬀective when l_danabasoglu =.true.

c. Anisotropic scheme

Behavior of GMANISOTROP option should be specified using the namelist listed on Table 9.13.
Table9.13 Namelist nml_gm_anisotrop

variable namecscl_isotrop

units
1

description
factor for anisotropy in GM diﬀusivity– 118 –

usage
if GMANISOTROPTECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 10

SGS parameterization of vertical mixing oftracers
This chapter explains subgrid-scale parameterizations of vertical mixing of tracers.

10.1 Vertical diﬀusion
Vertical diﬀusion term takes the form of Laplacian and the vertical diﬀusion flux is proportional to the vertical gradient oftracer. The finite diﬀerence form is as follows:
FZDi, j,k = −κ z (areat)i, j,k+ 1 δz Ti, j,k ,(10.1)

2

where the use of (areat)i, j,k+ 1 implies that the flux occurs only through the oceanic part of the grid interface and2δz Ti, j,k ≡Ti, j,k− 1 − Ti, j,k+ 1
2

∆zk

2

.

(10.2)

Note that, for simplicity, the change of the grid thickness at the bottom and fluctuations of the surface height are notconsidered in the grid distance ∆zk when calculating the gradient.
In most realistic simulations, a backward (implicit) scheme is used in the time integration (VVDIMP option; Section 19.5)
because high diﬀusivity is expected owing to the choice of parameterizations needed for realistic simulations. Otherwise,a forward scheme is used.

10.1.1 Specification of coeﬃcient
Background vertical diﬀusivity, which is horizontally uniform, a function of depth, and fixed in time, should be always
given. Additionally, a three dimensional distribution can be set by selecting VMBG3D option to incorporate locally enhanced
mixing processes induced by interaction between the bottom topography and tidal currents (e.g., St. Laurent et al., 2002).
With this choice, three dimensional distributions for vertical diﬀusivity and viscosity should be prepared in advance.
Tables 10.1 and 10.2 summarizes how to give background vertical diﬀusivity.
Table10.1 Namelist nml_diff_vert_bg. Specify only one of the two variables

variable name
diff_vert_bg_cm2psunits
cm2 s−1

file_diff_vert_1d_cm2ps

description
vertically uniform value of background vertical diﬀusivity
file having vertical 1D distribution

usage
Usable only if a vertically uniform
value is intended
cannot be specified with the above

Table10.2 Namelist nml_vmbg3d. Specify when VMBG3D is selectedvariable namefile_vmix_3d
imvm
jmvm

units

description
file having 3D distribution
east-west data size
north-south data size

usage

Continued on next page

– 119 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201710.2 Convective adjustment

variable namel_vmintpol
slatvmslonvmdlatvmdlonvmunits
logicalTable 10.2 – continued from previous page
description
usage
interpolate input data to model grid points
latitude of the western end
only if l_vmintpol = .true.
longitude of the southern end
only if l_vmintpol = .true.
uniform grid spacing in the meridional direction
only if l_vmintpol = .true.
uniform grid spacing in the zonal direction
only if l_vmintpol = .true.

In addition to the static background profiles, the following processes give varying vertical diﬀusivity coeﬃcients atevery model time step.
• Surface mixed layer models (TURBULENCE option).
• Vertical component of isopycnal diﬀusion (ISOPYCNAL option).
• Enhanced diﬀusivity (= 1.0 m2 s−1 = 104 cm2 s−1 ) where the stratification is unstable (DIFAJS option).
• Enhanced diﬀusivity around rivermouths to avoid too low salinity if river run-oﬀ is received by the model (RUNOFF
option). This scheme is especially needed when positive definiteness is not guaranteed by a tracer advection
algorithm. See Table 10.3 for how to specify the mixing.
The vertical diﬀusion for "this" time step is taken as the largest of the above estimations.
Table10.3 Namelist nml_vmix_river. Specify when RUNOFF is selectedvariable namel_enhance_vmix_rivmouth

units
logicaldiff_max_vmixriv_cm2ps

cm2 sec−1

depth_max_vmixriv_cm

cm

para_vmixriv_1

1

para_vmixriv_2

1

description
diﬀusivity is enlarged around therivermouth
maximum value of the enlarged diffusivity (= κ rivmax )vertical diﬀusion is enlarged fromsurface to this depth, this is also usedby subroutine salinity_limit
parameter for the enlarged vertical
diﬀusion formula (= a). See below.
parameter for the enlarged vertical diﬀusion formula (= b). Enlarged diﬀusion is calculated asκ riv = min(alog10 Wriver +b, κ rivmax ),where Wriver is river discharge ratein cm sec−1

usage
default = .false.
default = 1 × 104 cm2 sec−1
default = 30 × 102 cm
default = 10
default = 7

10.2 Convective adjustment
Convective adjustment is realized by replacing the density (temperature and salinity) that is statically unstable (the upper
density exceeds the lower density) in a water column with the averaged density between neighboring levels (neutralization),
considering that interior convection occurs in that place. Most of the realistic phenomena represented by the convective
adjustment include the developing mixed layer due to surface cooling during winter. Convective adjustment also includes
the case in which dense bottom water flows out the sill and flows down along the slope. Moreover, the convective adjustment
includes the practical eﬀect that it suppresses disturbances caused by the numerical calculation error and smooths the
distribution.
In general, there are three numerical schemes for convective adjustment. In the simplest one, adjustment is done for a pairof two neighboring levels, and then for a pair of another two neighboring levels. By repeating this procedure, it attempts
to neutralize the density in the unstable part. This procedure is simple at each step, but it has a defect that the finite-time
repetition does not necessarily guarantee reaching the complete averaged value. In the second scheme, adjustment isdone by assigning a high vertical diﬀusivity between the two levels that are statically unstable and by solving the vertical
diﬀusion term using an implicit method. This method cannot remove the unstable condition completely in one procedure.
However, it has good calculation eﬃciency for the case where the model has a high vertical diﬀusivity already due to themixed layer or isopycnal diﬀusion schemes and thus needs an implicit method to solve it. In MRI.COM, this scheme isinvoked by specifying DIFAJS option. The vertical diﬀusivity between the unstable grid points is set to 10 000 cm2 s−1 .
– 120 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 10

SGS parameterization of vertical mixing of tracers

In the third scheme, the unstable part is first neutralized. The stability at the top and bottom of the neutralized column isthen reexamined. If the unstable condition remains, the part including the already-neutralized column is re-neutralized.
This procedure continues until the instability at the top and bottom of the neutralized column disappears. This method
can remove the unstable part completely and thus is called "Complete Convection," but it requires a number of iterations,the vertical level size minus one, at maximum. The third method, which is the default scheme in MRI.COM, is explained
below (Ishizaki, 1997).

10.2.1 Algorithm
In order to minimize the judgment process ("IF" statement) and replace it by arithmetic calculation, this scheme defines
two integer indices, α k and βk , at the layer boundaries, and six real variables TU k , TLk , SU k , SLk , VU k , and VLk , (k =1, KM − 1), in addition to the vertical rows of temperature, salinity, and density Tk , Sk , Rk , (k = 12 , KM − 12 ) (KM is thenumber of levels; see Figure 10.1). The level at the vertical boundary of a T-cell corresponds to the integer k. The index
α k indicates an unstable part within a water column: α k = 1 if it is unstable at the level between k − 12 and k + 12 , andα k = 0 if it is neutral or stable. The index βk memorizes the mixed part: βk = 1 at the boundary where it is neutral asa result of mixing, and βk = 0 elsewhere. Variables TU k , SU k , and VU k and TLk , SLk , and VLk are temperature, salinityand volume accumulated by multiplying α above the level k and below the level k, respectively, and are expressed by thefollowing recursive relation.
VU 1 = α1V 1 ,2

VU 2 = α2 (V1+ 1 + α1 V 1 ) = α2 (V1+ 1 + VU 1 ),
2

2

2

· ··,
VU k = α k (Vk− 1 + VU k−1 ),2· ··,VU KM−1 = αKM−1 (VKM−1− 1 + VU KM−2 ),(10.3)2andVLKM−1 = αKM−1VKM− 1 ,2

VLKM−2 = αKM−2 (VKM−1− 1 + αKM−1VKM− 1 ) = αKM−2 (VKM−1− 1 + VLKM−1 ),
2

2

2

· ··,
VLk = α k (Vk+ 1 + VLk+1 ),2· ··,VL1 = α1 (V1+ 1 + VL2 ),(10.4)

2

where Vk+ 1 denotes a volume of the cell at the level k + 12 . In a similar way, other quantities are expressed as follows:2

TU 1 = α1T 1 V 1 ,2

TU k = α k (Tk− 1 Vk− 1 + TU k−1 ),SU 1 = α1 S 1 V 1 ,SU k = α k (Sk− 1 Vk− 1 + SU k−1 ),

2

2

2

2

2

2

2

TLKM−1 = αKM−1TKM− 1 VKM− 1 ,2

T L k = α k (Tk+ 1 Vk+ 1 + TLk+1 ),SLKM−1 = αKM−1 SKM− 1 VKM− 1 ,SL k = α k (Sk+ 1 Vk+ 1 + SLk+1 ),

2

2

2

2

2

2

where Tk+ 1 and Sk+ 1 are temperature and salinity at the level k + 12 .
22
According to this definition, if α k = 1 and elsewhere 0, we get
VU k + VLk = Vk− 1 + Vk+ 1 ,2

2

TU k + TLk = Tk− 1 Vk− 1 + Tk+ 1 Vk+ 1 ,
2

2

2

2

SU k + SLk = Sk− 1 Vk− 1 + Sk+ 1 Vk+ 1 ,2

2– 121 –

2

2

2

(10.5)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201710.2 Convective adjustment
Surface
α1

k=1/2

T1/2

k=1+1/2

T1+1/2

k-1/2

Tk-1/2

k+1/2

Tk+1/2

k=1

α2

k=2

αk-1

k-1

αk

k

αk+1

k+1

αKM-2
k=KM-2

TKM-1-1/2
k=KM-1-1/2

αKM-1
k=KM-1

TKM-1/2
k=KM-1/2
k=KM

BottomFigure10.1 Reference vertical grid points in Section 10.2

indicating a volume and accumulated temperature and salinity in an unstable part and
TU k + TLk,
VU k + VLk
SU k + SLk
SM k− 1 ,k+ 1 =,2
2
VU k + VLk

TM k− 1 ,k+ 1 =2

2(10.6)

are volume averaged temperature and salinity, respectively.
If the level k constitutes a series of the unstable part, the same equation holds for the averaged temperature and salinity.
For example, let α k−1 = α k = 1 and α k−2 = α k+1 = 0,
VU k−1 + VLk−1 = VU k + VLk= Vk−1− 1 + Vk− 1 + Vk+ 1 ,2

2

2

TU k−1 + TLk−1 = TU k + TLk= Tk−1− 1 Vk−1− 1 + Tk− 1 Vk− 1 + Tk+ 1 Vk+ 1 ,
2

2

2

2

2

2

SU k−1 + SLk−1 = SU k + SLk= Sk−1− 1 Vk−1− 1 + Sk− 1 Vk− 1 + Sk+ 1 Vk+ 1 ,
2

2

2

2

2

2

(10.7)

andTU k−1 + TLk−1
TU k + TLk=,
VU k−1 + VLk−1 VU k + VLk
SU k−1 + SLk−1
SU k + SLk
SM k−1− 1 ,k+ 1 ==.
2
2
VU k−1 + VLk−1 VU k + VLk

TM k−1− 1 ,k+ 1 =2

2

These are averages of the three layer, k − 1 − 12 , k − 12 , and k + 21 .
– 122 –

(10.8)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 10

SGS parameterization of vertical mixing of tracers

10.2.2 Numerical procedure
In summary, numerical procedures are summarized as follows.
[1] Density is calculated at the intermediate depth between adjacent levels using a pair of temperature and salinity andis judged to be statically stable or unstable. If an instability occurs, α(α 1 ) is replaced by 1, otherwise by 0. At this stage,β( β 0 ) is set to 0, where the superscript denotes the number of the iteration.
After this preprocessing, the following procedure (represented by n-th) is repeated until the instability is removed.
[2] Based on equations (10.3) to (10.5), VU, TU, SU, VL, TL, and SL are calculated using α n for a water column that
includes an unstable part.
[3] The vertical mean TM and SM are calculated for the unstable part using equation (10.6) and substituted for theoriginal values of T and S. This change modifies the density at the intermediate depth in [1].
[4] The value of α n is stored in β n . β n = 1 is set if α n = 1, or α n = 0 and β n−1 = 1, and otherwise β n = 0. This ispresented by the following:(10.9)
βkn = α kn + βkn−1 (1 − α kn ).
[5] The static stability is judged only for βkn = 0. Let α kn+1 = 1 if statically unstable, and 0 otherwise. If there is no
unstable part, the procedure for that water column is completed.
[6] For a water column which still includes an unstable part, modification for α kn+1 is done using βkn by the following.
After the procedure [2], any instability will be found only at the bottom of the part that is neutral as a result of prior mixing.
In that case, the neutral part must be treated as an unstable part, that is, α kn+1 = 1. On the other hand, no more procedureis needed if the upper and lower end is stable, giving α kn+1 = 0. This is done by a recursive formula going down and up in
the following.
γ1 = α1(n+1),

γk = α k(n+1) + (1 − α k(n+1) ) βk(n) γk−1

(n+1)
αKM−1
= γKM−1,(n+1)
α k(n+1) = γk + (1 − γk ) βk(n) α k+1,(10.10)where γ is a work variable, but may be treated as α itself in a FORTRAN program. Then, the procedure goes back to [2].
Table 10.4 shows an example of the case with six levels. Static instability is removed after the three-time iteration. The
second column of α in the table is the result of the corrected α kn+1 using βkn based on equation (10.10), as described in [6].
Note that βk0 = 0, though there is no description in the table.

– 123 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201710.2 Convective adjustment

Table10.4 Example of the convective adjustment procedure
n

1

VU+VL

TU+TL

β

V 1 + V1 1 + V2 12
2
2
V 1 + V1 1 + V2 1

T 1 V 1 + T1 1 V1 1 + T2 1 V2 1
2
2
2
2
2
2
T 1 V 1 + T1 1 V1 1 + T2 1 V2 1

1
1

0
0
V5 1

0
0
V4 1 + V5 1

00
T4 1 V4 1 + T5 1 V5 1

0
0
1

V1 1 + V2 1 + V3 1

V 1 + V1 1 + V2 1 + V3 1

VU

VL

1
2

1
1

1
1

V1
2
V 1 + V1 1

3
4
5

0
0
1

0
0
1

0
0
V4 1

1

0

1

V1

2

2

α

k

3

0
1

1
1

2

2

V1 1 + V2 1
2
2
V2 1
2

2

2

2

2

V 1 + V1 1
2

2

2

2

2

V2 1 + V3 12

V 1 + V1 1 + V2 1
2

2

2

2

V3 1

2

2

2

2

2

2

2

2

2

2

2

2

2

2

3
∑

2

2

2

2

2

Tk+ 1 Vk+ 1

1

Tk+ 1 Vk+ 1

1

Tk+ 1 Vk+ 1

1

0
0

0
1

2

k=0
3
∑

2

2

2

k=0
3
∑

2

V 1 + V1 1 + V2 1 + V3 12

2

2

V 1 + V1 1 + V2 1 + V3 1
2

2

2

2

2

2

k=0

4
5

0
0

0
0

0
0

0
0

1

0

1

V1

V1 1 + V2 1 + V3 1 + V4 1 + V5 1

2

3

3
4
5

0
0
1
0

1
1
1
1

2

2

V 1 + V1 1
2

2

2

2

2

2

2

2

2

2

2

2

2

2

V 1 + V1 1 + V2 1 + V3 1 + V4 1
2

2

2

2

2

Vk+ 1
Vk+ 1
Vk+ 1

2

Vk+ 1

Tk+ 1 Vk+ 1

1

Tk+ 1 Vk+ 1

1

Tk+ 1 Vk+ 1

1

Tk+ 1 Vk+ 1

1

2

k=0
5
∑

2

k=0
5
∑

2

k=0

1

2

k=0
5
∑

2

k=0
5
∑

Tk+ 1 Vk+ 1
2

k=0
5
∑

2

k=0
5
∑

2

– 124 –

5
∑

2

k=0
5
∑

2

V5 1

Vk+ 1
2

k=0
5
∑

2

V4 1 + V5 1
2

5
∑

2

V3 1 + V4 1 + V5 1

V 1 + V1 1 + V2 1 + V3 1
2

2

V2 1 + V3 1 + V4 1 + V5 1

V 1 + V1 1 + V2 1
2

0
0

2

k=0

2

2

2

2

2

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 11

Biogeochemical model
There are several options for Biogeochemical models in MRI.COM. These biogeochemical models have been developed forboth ocean-only and coupled ocean-atmosphere-vegetation carbon cycle studies. They feature an explicit representationof a marine ecosystem, which is assumed to be limited by light, temperature, and nutrients availability. This chapter
describes the details of the biogeochemical models.

11.1 Inorganic carbon cycle and biological model
Biogeochemical models are composed of inorganic carbon-cycle and ecosystem component models. In the inorganiccarbon-cycle component, the partial pressure of CO2 at the sea surface (pCO2 ) is diagnosed from the values of dissolved
inorganic carbon (DIC) and Alkalinity (Alk) at the sea surface, which should be aﬀected by the ecosystem component.
The diﬀerence in pCO2 between the atmosphere and ocean determines uptake or release of CO2 from the ocean to theatmosphere and is essential for simulating the CO2 concentration in the atmosphere. Inorganic carbonate chemistry
and partial pressure physics are well understood and can be reproduced with fair accuracy. The ecosystem component
deals with various biological activities, and gives sources and sinks of the nutrients, DIC, Alk, and dissolved oxygen
through these activities. Our knowledge of these activities is far from complete, and they are diﬃcult to estimate even instate-of-the-art models.
There are many biological models and methods for calculating the ecosystem components. One of the simplest biologicalmodels has only one nutrient component (such as PO4 ) as a prognostic variable and calculates neither phytoplankton nor
zooplankton explicitly. In these cases, the export of biologically generated soft tissue (organic matter) and hard tissue(carbonate) to the deep ocean, collectively known as the biological pump, is parameterized in terms of temperature,salinity, shortwave radiation, and nutrients.
A Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) model is more complex than the above model, but still a
simple biological model. The NPZD model has four prognostic variables (nutrient, phytoplankton, zooplankton, anddetritus). Though parameterized in a simple form, basic biological activities, such as photosynthesis, excretion, grazing,and mortality are explicitly calculated.
More complex models classify phytoplankton and zooplankton into several groups, and deal with many complex
interactions between them. In general, it is expected that the more complex the biological model becomes the more
realistic pattern the model can simulate. However, because of our incomplete knowledge about the biological activities,
the complex models do not always yield better results, even though they require more computer resources.
To simulate the carbon cycle in the ocean, some biological processes should be calculated in the ecosystem component
to obtain DIC at the sea surface. However, the carbon cycle component is not always necessary when our interests are tosimulate the ecosystem itself. The Ocean Carbon-Cycle Model Intercomparison Project (OCMIP) protocols and studies
of Yamanaka and Tajika (1996) and Obata and Kitamura (2003) focus on the former carbon cycle in the ocean, and theecosystem components in these studies are quite simple. Biogeochemical models adopted in MRI.COM are classified inthis category. The latter studies usually use complex biological models such as NEMURO (Kishi et al., 2001). Of course,
this type of model could be adopted as an ecosystem component of the biogeochemical model in the former studies inhopes of better simulation of carbon cycle.
Originally, the carbon cycle component followed the OCMIP protocols (Orr et al., 1999) whose authority is recognizedin the community. Recently, biogeochemical protocols for CMIP6 require that the former OCMIP code should be replacedby mocsy routines (Orr and Epitalon, 2015) to use the equilibrium constants recommended for best practices (Orr et al.,2016). MRI.COM can choose either option by setting a namelist. Here we describe the procedure using the mocsy routines.
MRI.COM has several options for the ecosystem component. At present, MRI.COM can incorporate the Obata andKitamura model (Obata and Kitamura, 2003) or an NPZD model based on Oschlies (2001). The biogeochemical modelof MRI.COM is largely based on Schmittner et al. (2008) when an NPZD model is adopted as an ecosystem component.
Basic units in MRI.COM are cgs, but in these biogeochemical subroutines, we use MKS units for the sake of future– 125 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201711.2 Governing equations
development. We use mol/m3 for the units of nutrients. When the coeﬃcients of their model are applied, they should be
converted to the corresponding units.

11.2 Governing equations
Here we describe the biogeochemical models of MRI.COM. When an NPZD model is incorporated as the ecosystemcomponent, the governing equations are as follows. When Obata and Kitamura model is used instead of the NPZD model,
the first four biogeochemical compartments (DIC, Alk, PO4 , and O2 ) are used.
∂DIC= −A(DIC) + D (DIC) + Sb (DIC) + Jv (DIC) + Jg (DIC),(11.1)∂t∂Alk= −A(Alk) + D(Alk) + Sb (Alk) + Jv (Alk),(11.2)∂t∂[PO4 ]= −A([PO4 ]) + D([PO4 ]) + Sb ([PO4 ]),(11.3)∂t∂[O2 ]= −A([O2 ]) + D ([O2 ]) + Sb ([O2 ]) + Jg ([O2 ]),(11.4)∂t∂[NO3 ]= −A([NO3 ]) + D([NO3 ]) + Sb ([NO3 ]),(11.5)∂t∂[PhyPl]= −A([PhyPl]) + D([PhyPl]) + Sb ([PhyPl]),(11.6)∂t∂[ZooPl]= −A([ZooPl]) + D([ZooPl]) + Sb ([ZooPl]),(11.7)∂t∂[Detri](11.8)= −A([Detri]) + D ([Detri]) + Sb ([Detri]),∂twhere A() and D() represent advection and diﬀusion operator, respectively, and Sb () is source minus sink due to thebiogeochemical activities. The square brackets mean dissolved concentration in mol/m3 of the substance within them.
The terms represented by Jg () and Jv () are the air-sea gas fluxes (Section 11.3.1) and the dilution and concentration eﬀects
of evaporation and precipitation on DIC and Alk (Section 11.3.2), respectively, which appear only at the sea surface. The
term Jg () is calculated based on the OMIP protocol by using the air-sea gas transfer velocity and concentration in theseawater. The term Jv () appears only when the salinity flux is given virtually instead of the increase or decrease of thevolume at the surface layers due to evaporation and precipitation.

11.3 Carbon cycle component
To estimate Jg and Jv , we follow the protocol of OMIP, which is described in detail in Orr et al. (2016). The program tocalculate them is based on the mocsy subroutines. We have modified this subroutine so that it can be used in the calculationof the MRI.COM code.

11.3.1 Air-sea gas exchange fluxes at the sea surface ( Jg )The air-sea gas transfer must be calculated for DIC and [O2 ]. The terms Jg (DIC) and Jg ([O2 ]) appear only in theuppermost layer. When these fluxes are expressed as Fg (DIC) and Fg ([O2 ]), Jg (DIC) and Jg ([O2 ]) are given as follows:
Jg (DIC) =Fg (DIC),∆z 1(11.9)

2

Fg ([O2 ])Jg ([O2 ]) =,∆z 1

(11.10)2

where
Fg (DIC) = KwCO2 ∗ ([CO2 ]sat − [CO2 ]surf ),(11.11)Fg ([O2 ]) = KwO2 ∗ ([O2 ]sat − [O2 ]surf ),(11.12)and ∆z 1 is the thickness of the first layer of the model. Here a standard gas transfer formulation is adopted. KwCO2 and KwO2
2
are their gas transfer velocity, [CO2 ]sat and [O2 ]sat are their saturation concentrations with respect to the atmosphere, and[CO2 ]surf , [O2 ]surf , are their surface concentration.
– 126 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 11

Biogeochemical modela. Gas transfer velocity

The coeﬃcients KwCO2 and KwO2 are the air-sea gas transfer (piston) velocity and are diagnosed as follows:() −1/2ScCO2
2U10(1 − f i ),660
( O2 ) −1/2
Sc
O2
2
U10
(1 − f i ),Kw = a660

KwCO2 = a(11.13)(11.14)where• f i is the fraction of the sea surface covered with ice,
• U10 is 10 m scalar wind speed,• a is the coeﬃcient of 0.251 cm hr−1 /(m s−1 ) 2 , which will give Kw in units of cm hr−1 when winds speeds (U10 ) arein m s−1 . This is specified in the OMIP protocol,• ScCO2 and ScO2 are the Schmidt numbers for CO2 and O2 .
b. Saturation concentration with respect to the atmosphere

The surface water gas concentration in equilibrium with the atmosphere (saturation concentrations) for gas A has thefollowing relationship.
[A]sat = K0 f A= K0 C f p A= K0 C f (Pa − pH2 O)x A= ϕAx A
Pa
(with errors less than 0.1%)≃ ◦ ϕ0A x APaPa= ◦ [A]0sat,
Pa(11.15)(11.16)(11.17)(11.18)(11.19)(11.20)where K0 is its solubility, f A is its atmospheric fugacity, C f is its fugacity coeﬃcient, p A is its atmospheric partial pressure,x A is its mole fraction in dry air, Pa is total atmospheric pressure, Pa◦ is the standard atmosphere (=1013.25 hPa), pH2 O
is water vapor pressure at saturation, ϕ A is its solubility function, ϕ0A and [A]0sat are its solubility function and saturation
concentration at the reference pressure.
Orr et al. (2016) explain two approaches to calculate the saturation concentration.
The first approach uses equations (11.19) or (11.20), computing solubility function ϕ0A or saturation concentration [A]0sat
at the reference atmospheric pressure first, and then converting them to those at the atmospheric pressure Pa . For gasesthat are often used as tracers in oceanography such as CFCs, ϕ0A can be expressed as a function of in-situ temperature andsalinity by using empirical fit. For oxygen, an empirical fit for [O2 ]0sat is available. This is a rather conventional approach,and is used in MRI.COM except for carbon.
The other approach uses equations (11.15)-(11.17), computing its partial pressure p A ≡ (Pa −pH2 O)x A, then multiplying
by K ′ ≡ K0 C f , or multiplying by C f and then multiplying by K0 if atmospheric fugacity ( f A ) is needed. For typicalgases, K ′, K0 , and C f are available as a function of temperature and salinity. The mocsy routine adopts this approach, and
MRI.COM also follows it for calculating carbon flux at the surface.
c. Computing CO2 concentrations at the surface

In the mocsy routine, [CO2 ]surf is diagnosed every step from DIC, Alk, temperature, salinity, [PO4 ], and silicate concentration at the surface. After [CO2 ]surf is computed, fugacity ( f CO2 ) and partial pressure (pCO2 ) of the ocean surface arecomputed as
f CO2 = [CO2 ]surf /K0,
pCO2 = f CO2 /C f ,where K0 and C f are estimated by using temperature and salinity.
– 127 –

(11.21)(11.22)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201711.3 Carbon cycle component
Diagnosis of [CO2 ]surf is the most complex of the above calculations and has the largest computational cost. To be
precise, diagnosis of [CO2 ] actually means diagnosing [CO2 ] + [H2 CO3 ], which are diﬃcult to distinguish analytically.
These two species are usually combined and the sum is expressed as the concentration of a hypothetical species, [CO∗2 ] or
[H2 CO∗3 ]. Here, the former notation is used. The relationship between this [CO∗2 ] and DIC is as follows:
DIC = [CO2 ] + [H2 CO3 ] + [HCO−3 ] + [CO2−3 ]= [CO∗2 ](11.23)+ [HCO−3 ] + [CO2−
3 ].

(11.24)In the OMIP protocol, the following equations are solved to obtain [CO∗2 ].
The equilibrium constants for dissociation reactions are:
K1 =[H+ ][HCO−3 ][CO∗2 ]

[H+ ][H2 PO−4 ][H3 PO4 ]

K2P =

KSi =

[H+ ][CO32− ],[HCO−3 ](11.25)[H+ ][B(OH) −4 ],[B(OH) 3 ]

KB =

K1P =

K2 =

[H+ ][HPO2−4 ][H2 PO−4 ](11.26)

K3P =[H+ ][PO3−4 ][HPO2−
4 ][H+ ][SiO(OH) −3 ],[Si(OH) 4 ],(11.27)(11.28)KW = [H+ ][OH− ],(11.29)[H+ ]F [SO2−4 ],−[H SO4 ](11.30)[H+ ]F [F− ],[HF](11.31)

KS =
and

KF =

where [H+ ] is the hydrogen ion concentration in sea water and [H+ ]F is the free hydrogen ion concentration. There isanother scale for the hydrogen ion concentration, the total hydrogen ion concentration [H+ ]T . The subscript T means"total" and F means "free." These three hydrogen ion concentrations are related as follows:
(
)
ST
FT
+
+
[H ] = [H ]F 1 ++,(11.32)
KS KF
(
)
ST
and
[H+ ]T = [H+ ]F 1 +.
(11.33)KS
There are three pH scales corresponding to these three hydrogen ion concentrations.
The equilibrium constants K x are given as a function of temperature, salinity, and pH. Note that the equilibrium constants
are given in terms of concentrations, and that all constants are referenced to the seawater pH scale, except for KS , whichis referenced to the free pH scale.
The total dissolved inorganic carbon, boron, phosphate, silicate, sulfate, and fluoride are expressed as follows:
DIC = [CO∗2 ] + [HCO−3 ] + [CO2−3 ],(11.34)

BT = [B(OH) 3 ] + [B(OH) −4 ],(11.35)

3−
PT = [H3 PO4 ] + [H2 PO−4 ] + [HPO2−
4 ] + [PO4 ],(11.36)

Si T = [Si(OH) 4 ] + [Si(OH) −3 ],(11.37)– 128 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 11

Biogeochemical modelandST = [HSO−4 ] + [SO2−
4 ],(11.38)

FT = [HF] + [F− ].

(11.39)Alkalinity used in this calculation is defined as follows:
−
−
2−
3−
−
Alk = [HCO−3 ] + 2[CO2−3 ] + [B(OH) 4 ] + [OH ] + [HPO4 ] + 2[PO4 ] + [SiO(OH) 3 ]

− [H

+

(11.40)]F − [HSO−4 ] − [HF] − [H3 PO4 ].

These expressions exclude the contribution of NH3 , HS− , and S2− .
If we assume that DIC, Alk, [PT ], and [SiT ] are known, this system contains 18 equations with 18 unknowns, so they canbe solved. Detailed solver method has been described in Orr and Epitalon (2015). The concentration [SiT ] is not predictedin the biogeochemical model adopted in MRI.COM but is rather specified by the annual mean value from WOA2013.

11.3.2 Dilution and concentration eﬀects of evaporation and precipitation on DIC and Alk
The dilution and concentration eﬀects of evaporation and precipitation significantly impact the concentrations of somechemical species in seawater. This is particularly true for DIC and Alk, which have large background concentrations
compared with their spatial variability. MRI.COM uses a free surface, so the impact of evaporation and precipitation isstraightforward to model unless SFLUXW or SFLUXR option is used. In these options, salinity flux is diagnosed and appliedinstead of the freshwater flux. In this case, the dilution and concentration eﬀect of evaporation (E) and precipitation (P)should be taken into account. Here, they are parameterized as virtual DIC and Alk fluxes, similar to the virtual salt flux
used in physical ocean GCMs.
In MRI.COM, the tendency of salinity due to the virtual salt flux is given bysflux(i, j) = −(P − E) ∗ S(i, j, 1)/∆z,(11.41)

where S(i, j, 1) and ∆z are the salinity and thickness of the uppermost layer. Note that the variable sflux(i, j) is not thesalinity flux but the time change rate of salinity due to the flux even though the spelling brings up the image of the flux.
In MRI.COM, DIC and Alk are modified by the virtual salt flux as follows:
Jv (DIC(i, j, 1)) = sflux(i, j)/S(i, j, 1) ∗ DIC(i, j, 1),
Jv (Alk(i, j, 1)) = sflux(i, j)/S(i, j, 1) ∗ Alk(i, j, 1).

(11.42)(11.43)

Strictly speaking, air-sea fluxes of fresh water impact other species. However, these modifications are not usually applied
because their spatial variabilities are significantly greater than those of DIC and Alk.
In the OMIP protocol, as well as the previous OCMIP protocol, the global averaged salinity Sg is used instead ofS(i, j, 1) in equations(11.42,11.43). In addition, globally integrated Jv (DIC) and Jv (Alk) are set to 0. In MRI.COM, thesemodifications are not the default considering the use in regional ocean models.

11.4 Obata and Kitamura modelThis section was contributed by A. Obata.
The Obata and Kitamura model used in MRI.COM simply represents the source and sink terms of DIC, Alk, [PO4 ], and[O2 ] due to the biogeochemical activities: new production driven by insolation and phosphate concentration in the surface
ocean, its export to depth, and remineralization in the deep ocean. According to the Michaelis-Menten kinetics (Dugdale,1967), phosphorus in the new production exported to depth (ExprodP) is parameterized as r L[PO4 ]2 /([PO4 ] + k), wherer is a proportional factor (r = 0.9 yr−1 ), L is the insolation normalized by the annual mean insolation on the equator, andk is the half-saturation constant (k = 0.377 mmol/m−3 ). The values of r and k are adjusted to reproduce the optimum
atmospheric CO2 concentration and ocean biogeochemical distribution for the preindustrial state of the model. The
relationship between the changes in the chemical composition of seawater and the composition of particulate organic
matter (POM) is assumed to follow the Redfield ratio P : N : C : O2 = 1 : 16 : 106 : −138 ≡ 1 : Rnp : Rcp : −Rop ,where Rab = A/B and "o" represents O2 (Redfield et al., 1963). The rain ratio of calcite to particulate organic carbon(POC) is 0.09, which is in the range proposed by Yamanaka and Tajika (1996). The surface thickness where the export
production occurs is fixed at 60 m. The vertical distribution of POM and calcite vertical flux below a depth of 100 m is
– 129 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201711.5 NPZD model
proportional to (z/100m) −0.9 and exp(−z/3500m) (z is the depth in meters), respectively, following the work of Yamanaka
and Tajika (1996). The remineralization of POM (RemiP for phosphorus) and the dissolution of calcite (SolnCa) at depth
are parameterized by these fluxes. Oxygen saturation is prescribed at the sea surface. The solubility of oxygen is computedfrom the formula of Weiss (1970). Source and sink terms of Sb () representing the above processes are as follows:
Sb (DIC) = Rcp ∗ RemiP + SolnCa − Rcp ∗ ExprodP
Sb (Alk) = 2 ∗ SolnCa + Rnp ∗ ExprodP − Rnp ∗ RemiPSb ([PO4 ]) = RemiP − ExprodP
Sb ([O2 ]) = − Rop ∗ Sb ([PO4 ])(11.44)(11.45)(11.46)(11.47)

11.5 NPZD model
The NPZD model used in MRI.COM is constructed on the assumptions that the biological elemental composition ratio isnearly constant (Redfield ratio) and that the concentration of organisms can be estimated by nitrogen or phosphorus. The
prognostic variables of nitrate (NO3 ), phytoplankton (PhyPl), zooplankton (ZooPl), and detritus (Detri) are normalized interms of nitrogen 1 mol/m3 . For example, [PhyPl] represents the concentration of phytoplankton estimated in terms ofnitrogen in one cubic meter (N mol/m3 ). The increase and decrease of carbon can be diagnosed by multiplying by Rcn .
Source and sink terms Sb () calculated in the NPZD model are as follows. Those for DIC and Alk, Sb (DIC) and Sb (Alk),
used for calculating the carbon cycle, are described later in this section.
Sb ([PhyPl]) = Priprod − MortP1 − MortP2 − GrP2Z
Sb ([ZooPl]) = assim ∗ GrP2Z − Excrtn − MortZ
Sb ([Detri]) = [(1 − assim) ∗ GrP2Z + MortP2 + MortZ] − RemiD − wdetri
Sb ([NO3 ]) = MortP1 + Excrtn + RmeiD − PriProd
Sb ([PO4 ]) = Rpn ∗ Sb ([NO3 ])Sb ([O2 ]) = − Ron ∗ Rnp ∗ Sb ([PO4 ])(11.48)(11.49)∂Detri∂z(11.50)(11.51)(11.52)(11.53)

There is no input from the atmosphere such as nitrogen fixation in the above equations, so the sum of these five equations
becomes zero at each grid point except for the term for detritus sinking (−wdetri ∂Detri∂z ). The term for detritus sinking
expresses the biological pump, whose role is to remove nutrients from the upper layers and transport them into the deep
ocean where the plankton cannot use the nutrients. When vertically integrated, the sum of each grid is 0 even though
this sinking term is included. The nutrients are transported horizontally through physical processes such as advection anddiﬀusion.
In general, the nitrate limitation is more severe than the phosphate limitation so it is not always necessary to calculate
phosphate. However, in the simpler model of Obata and Kitamura (2003), phosphate is used as a prognostic variable.
So, to be consistent, phosphate is calculated in the ecosystem component of MRI.COM. Next, we elaborate on the above
equations.

11.5.1 Description of each term
• Priprod = J (I, N, P) ∗ [PhyPl]
Primary production expresses photosynthesis (described in detail in the next subsection).
• MortP1 = ϕ P ∗ [PhyPl]
The conversion of mortality phytoplankton directly into nutrients. This term was introduced by Oschlies (2001) toincrease the primary production of subtropical gyre, where the nutrient limit is severe.
• MortP2 = ϕPP ∗ [PhyPl]2
The conversion from phytoplankton to detritus (normal mortality of phytoplankton).
• GrP2Z = G(P) ∗ [ZooPl]2
The grazing of zooplankton. There are a number of parameterizations of grazing. In this formulation, this is given
as G(P) = g ∗ ϵ ∗ [PhyPl]2 /(g + ϵ ∗ [PhyPl]2 ). Among the grazing, the ratio assim is used for the growth ofzooplankton, and the remainder (1 − assim) is converted to detritus.
• Excrtn = d ∗ [ZooPl]
Excretion of zooplankton. The excretion is dissolute and directly returned to nutrients (NO3 ).
– 130 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 11

Biogeochemical model• MortZ = ϕ Z ∗ [ZooPl]
The conversion from the zooplankton to detritus (mortality of zooplankton).
• RemiD = ϕ D ∗ [Detri]Remineralization of detritus. This is converted to nutrients through the activity of bacteria.
2TKOCT[2TQFWEVKQP

㧺

2
[

OQTVCNKV[

P

OQ

TV

KQ
GV

CN

ET

KV

GZ

TGOKPGTCNK\CVKQP&

ITC\KPI<GZRQTVOQTVCNKV[

Figure11.1 Schematic of NPZD model

11.5.2 Primary Production
The growth rate of phytoplankton is limited by the irradiance (I) and nutrients. This limitation is expressed in several
ways. Here we adopt an expression with a minimum function:
J (I, N, P) = min (JI , JN , JP ) ,(11.54)

where JI denotes the purely light-limited growth rate, and JN and JP are nutrient-limited growth rates that are functionsof nitrate or phosphate.
The light-limited growth is calculated as follows:
JI = [

Jmax αI]

2 + (αI) 2 1/2
Jmax.

(11.55)Here, Jmax is the light-saturated growth, which depends on temperature based on Eppley (1972) as
Jmax = a · bcθ ,(11.56)

where a = 0.6 day−1 , b = 1.066, and c = 1 (◦C) −1 . Note that the default values in MRI.COM are based on Schmittner
et al. (2008) and diﬀer from these values (see Table 11.1). Equation (11.55) is called Smith-type growth. The coeﬃcient
α in the equation is "the initial slope of photosynthesis versus irradiance (P-I) curve," that is,∂ JI.
I→0 ∂Iα = lim(11.57)

Thus, it represents how sensitive JI is to the irradiance when the light is weak.
Irradiance (I) depends on the angle of incidence and the refraction and absorption in the seawater.
()∫ z̃
I = Iz=0 PAR exp −k w z̃ − k ePdz ,0– 131 –

(11.58)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201711.5 NPZD model
where Iz=0 denotes the downward√shortwave radiation at the sea surface, PAR is the photosynthetically active radiation
ratio (0.43) and z̃ = z/ cos θ = z/ 1 − sin2 θ/1.332 is the eﬀective vertical coordinate (positive downward) with 1.33 as
the refraction index according to Snell’s law relating the zenith angle of incidence in air (θ) to the angle of incidence inwater. The angle of incidence θ is a function of the latitude ϕ and declination δ.
For the nutrient-limited growth rate (JN and JP ), we adopt the Optimal Uptake (OU) equation instead of the classic
Michaelis-Menten (MM) equation. For the classic MM equation, the nitrate-limited growth rate is expressed as
JN = JMM =Jmax N,
KN + N

(11.59)where K N is a half-saturation constant for NO3 uptake. In contrast, the Optimal Uptake (OU) equation for a nitrate isexpressed as follows:
V0 N
JN = JOU =,(11.60)√N + 2 VA00 N + VA00
where A0 and V0 are the potential maximum values of aﬃnity and uptake rate, respectively (see Smith et al. (2009) for
details). Optimal Uptake (OU) kinetics assumes a physiological trade-oﬀ between the eﬃciency of nutrient encounter at
the cell surface and the maximum rate at which a nutrient can be assimilated (Smith et al., 2009). The key idea is thatphytoplankton alters the number of its surface uptake sites (or ion channels), which determines the encounter timescale,versus internal enzymes, which assimilate the nutrients once encountered.
We set parameters V0 and A0 so that the rates of uptake, JMM and JOU , are equal at N = K N . In addition, we fix theratio V0 /A0 = αOU , where αOU is determined from fitting the data. This requires
(
)2
√
αOU
V0 = 0.5 1 +Jmax .
KN(11.61)

V0 N.
N + 2 αOU N + αOU(11.62)

Finally, we obtain
JOU =

√

We use αOU = 0.19, which is determined from the fitting of log K N vs log N in the wide range of N by Smith et al. (2009).

11.5.3 Variation of DIC and Alk due to biological activity
Production of DIC and Alk is controlled by changes in inorganic nutrients and calcium carbonate (CaCO3 ), in molar
numbers according to
Sb (DIC) = Sb ([PO4 ])Rcp − Sb ([CaCO3 ]),(11.63)
Sb (Alk) = −Sb ([NO3 ]) − 2 · Sb ([CaCO3 ]).

(11.64)Thus, only these source and sink terms of DIC and Alk are estimated. Since [PO4 ] and [NO3 ] are prognostic variables,
their source and sink are explicitly calculated by the biological model. In contrast, the downward movement of CaCO3 ismuch faster than the modeled downward velocity of water mass, so [CaCO3 ] is not a prognostic variable, and its source(Pr) and sink (Di) are diagnosed by the following equation,
Sb ([CaCO3 ]) = Pr ([CaCO3 ]) − Di([CaCO3 ]).

(11.65)

Following Schmittner et al. (2008), the source term (Pr ([CaCO3 ])) of calcium carbonate is determined by the productionof detritus as follows:
Pr ([CaCO3 ]) = [(1 − assim) ∗ [GrP2Z] + [MortP2] + [MortZ]] RCaCO3 /POC RC:N ,(11.66)

where assim, GrP2Z, MortP2, and MortZ are as described above. The sink term (Di([CaCO3 ])) of calcium carbonate isparameterized as∫)d (exp (−z/DCaCO3 ) ,(11.67)
Di([CaCO3 ]) =Pr ([CaCO3 ])dz ·dzwhich expresses an instantaneous sinking with an e-holding depth of DCaCO3 = 3500 m. In this equation, z is positive
downward. This depth of 3500 m was estimated by Yamanaka and Tajika (1996) to reproduce the observed nutrient profile.
– 132 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 11

Biogeochemical modelThis value is standard and also is used in the simple biological model in the protocol of OCMIP. The vertical integral of
the source minus sink should be zero. Thus, when the sea bottom appears before the sum becomes zero, the remaining
calcium carbonate is assumed to be dissolved in the lowermost layer. By using the ratio RCaCO3 /POC = 0.035 used by
Schmittner et al. (2008), the resultant global mean Rain ratio should be roughly consistent with the recently estimatedrange (0.07 to 0.11) based on various observations.

11.6 Usage
In MRI.COM, carbon, dissolved oxygen, and ecosystem model as well as other passive tracers like CFCs can be calculated
together or separately. Options described in configure.in determines what components should be used in MRI.COM.
The combination of the options also determines the number of passive tracers that is necessary for calculation. See detailfor the program tracer_vars.F90.
• When CARBON is used, carbon component is calculated in the model. For the ocean only model, the CBNHSTRUN
option should be used to set the atmospheric CO2. Its internal source or sink is calculated in an ecosystem model.
• When CARBON is used, dissolved oxygen component is calculated.
• When CBNHSTRUN is used, the ‘AtmosphericxCO2’ in ppm should be applied as additional atmospheric forcing.
This is the standard configuration for ocean only run. See Chapter 14 for the detailed setting.
• When O2 is used, dissolved oxygen component is calculated. Its internal source or sink is calculated in an ecosystemmodel.
For ecosystem model, MRI.COM can use NPZD or OBT options.
• When NPZD is used, an NPZD model is used as the ecosystem component. In addition, when CHLMA94 option
is used, the chlorophyll concentration is considered to calculate the shortwave penetration following Morel andAntoine (1994).
The parameters of NPZD are set in namelist nml_bioNPZD. The default values are based on Schmittner et al. (2008)
and listed on Table 11.1. If the parameters of Oschlies (2001) are used, the high nutrient-low chlorophyll (HNLC)
region in the North Pacific is not appropriately expressed. This may be because the parameters of Oschlies (2001)
are calibrated for the North Atlantic biological model. The commonly used unit of time in biological models is[day]. Thus, in the namelist, the time unit of the biological parameter is specified by using the unit [day]. In themodel, the time unit is converted to seconds, [sec].
• When OBT is used, a simple biological model of Obata and Kitamura (2003) is applied as the ecosystem component.
Restart files for the passive tracers should be specified by repeatedly writing nmlrs_ptrc. The order of passive tracers
is determined in tracer_vars.F90. Following is an example when CARBON, O2, and NPZD options are chosen. In thiscase, numtrc_p = 8 should be specified in configure.in.
An example specification of restart files when CARBON, O2, and NPZD options are chosen.
&nmlrs_ptrc fname=’result/rs_dic’/&nmlrs_ptrc fname=’result/rs_alk’/&nmlrs_ptrc fname=’result/rs_o2’ /&nmlrs_ptrc fname=’result/rs_po4’ /&nmlrs_ptrc fname=’result/rs_no3’ /&nmlrs_ptrc fname=’result/rs_PhyPl’/&nmlrs_ptrc fname=’result/rs_ZooPl’/&nmlrs_ptrc fname=’result/rs_Detri’/– 133 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201711.7 Program structure
Table11.1 Parameters used for the NPZD ecosystem component (NPZD).

variable namealphabio
abio
bbio
cbio
PARbiodkcbiodkwbiork1bioNO3
rk1bioPO4alpha_ougbio
epsbiophiphyphiphyqa_npz
phizood_npz
reminaw_detr
fac_wdetr

c_mrtnRcn
Ron
Rnp
dp_euphdp_mrtndp_eprdc
Rcaco3poc
Dcaco3
shwv_intv

description

units

Initial slope of P-I curve
Maximum growth rate parameter
Maximum growth rate = a * b ** (c * T)
Photosynthetically active radiation
Light attenuation due to phytoplankton
Light attenuation in the water
Half-saturation constant for NO3 uptake
Half-saturation constant for PO4 uptake
Fitting constant for Optical Uptake kinetics
Maximum grazing rate
Prey capture rate
Specific mortality/recycling rate
Quadratic mortality rate
Assimilation eﬃciency
Quadratic mortality of zooplankton
Excretion
Remineralization rate
Sinking velocity
Arbitrary parameter for numerical stability.
When the concentration of detritus in the n+1 st level is higher
than fac_wdetr times that in the n th level, w_detr is set to
0 between n and n+1 level.
Dimensionless scaling factor for Martin et al. (1987)
Phi(z) = Phi(zo) * (z/dp_mrtn)**(-c_mrtn)
Molar elemental ratio (C/N)
Molar elemental ratio (O2 /N)
Molar elemental ratio (N/P)
Maximum depth of euphotic zone
Characteristic depth of martin curve
The depth where the bio-export is diagnosed.
This value should be less than dp_mrtn.
CaCO3 over nonphotosynthetic POC production ratio
CaCO3 remineralization e-folding depthInterval for calculating the irradiance and light-limited growthrate. This must be a divisor of the time step for tracer.

11.7 Program structure
ogcm__ini|+-- ptrc_ctl__ini|||+-- ptrc__ini|||+-- cbn__ini|||+-- cbn__ini_history|+-- tracer__ctl__ini|+-- tracer__ini|+-- bio__ini– 134 –

(W m−2 ) −1 day−1
day−1

m−1 (mol m−3 ) −1
m−1
mol m−3
mol m−3
day−1(mol m−3 ) −2 day−1s−1−3(mol m ) −1 day−1(mol m−3 ) −1 day−1
day−1
day−1
m day−1

default value
0.1d0
0.2d0
1.066d0
1.d0
0.43d00.03d30.04d00.7d-30.0d0
0.19d01.575d0
1.6d6
0.014d00.05d30.925d0
0.34d30.01d00.048d0
2.0d0
3.d0

0.858d0

m
m
m

7.d0
10.d0
16.d0
150.d0400.d0126.d0m
min

0.05d03500.d0
10.d0

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

Biogeochemical model|+-- bio__ini_history

ogcm__run|+-- part_1|||+-- ptrc_ctl__main
|
|
|
|
|
+-- ptrc__surfflux
|
|
|
|
|
|
|
+-- cbn__set_xco2a
|
|
|
|
|
|
|
+-- cbn__calc_virtual_flux
|
|
|
|
|
|
|
+-- cbn__calc_co2_flux
|
|
|
|
|
|
|
+-- co2flux_mocsy
|
|
|
|
|
+-- ptrc__predict_surf
|
|
|
|
|
+-- cbn__predict_surf|||+-- tracer_ctl__predict|||+-- tracer__internal_source|||+-- bio__predict|+-- hist_ctl|+-- hist_ctl__write|+-- cbn__write_history– 135 –

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 12

Inert tracersThis chapter explains inert tracers implemented in MRI.COM. In general, to include an inert tracer in a model integration,
you need to specify the model options and tracer attributes corresponding to that tracer. General explanation on how to
specify attributes of a tracer is given at Section 13.3.

12.1 Ideal age tracer
Ideal age of a water mass is the time in year since it last contacted with the sea surface. This tracer is introduced byEngland (1995).

12.1.1 Source / Sink term
Source and sink of an age tracer a [year], is expressed as follows. In the oceanic interior,
1∂a= −A(a) + D(a) +∂t
the total number of seconds of this year

year/sec,(12.1)

where A() and D() represent advection and diﬀusion operator, respectively. Note that the r. h. s. depends on whether it
is leap year or not. The tracer ages slowly in leap years than in normal years.
At the sea surface, the ideal age is set to zero:a(k = 1) = 0.

(12.2)

12.1.2 Usage
Required Model option IDEALAGE
Required Add one to numtrc_p
Required Restart file for passive tracer (nmlrs_ptrc) (See Section 11.6)
Required Name for nml_tracer_data is ’Ideal Age Tracer’ (See Section 13.3)
Optional Namelist nml_tracer_idealage_start to specify the base date and time [year, month, day, hour, minute,second] for the age of water.

12.2 CFCs
MRI.COM can treat CFCs (CFC11 and CFC12) following the protocols of OMIP (Orr et al., 2016). There is no sourceand sink for CFCs in the interior.

12.2.1 Surface boundary condition
The CFCs have air-sea gas fluxes at the sea surface as the source. Their dissolved concentrations [A] [mol m−3 ] evolve
according to the advection (A()) and diﬀusion (D()) in the oceanic interior,∂[A]= −A([A]) + D([A]) + Jg ([A]).
∂t– 137 –

(12.3)TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201712.3 SF6
Jg () represents source/sink due to the air-sea gas fluxes F ([A]) at the sea surface,
Jg ([A]) =F ([A]),∆z 1(12.4)

2

withF ( A) = k w ([A]sat − [A]),(12.5)

where k w is its gas transfer velocity, [A]sat is the surface gas concentration in equilibrium with the atmosphere at an
atmospheric pressure at the surface (Pa ). k w is a function of wind velocity, surface temperature, and atmospheric pressure
at the sea surface. [A]sat is given by[A]sat = ϕ0A x A,
(12.6)
where x A is its mole fraction in dry air. The combined solubility term ϕ0A is computed using the empirical fit of temperatureand salinity. See Orr et al. (2016) for further details.

12.2.2 Usage
Required Model option CFC
Required Add two to numtrc_p
Required Restart file for passive tracer (nmlrs_ptrc) (See Section 11.6)
Required Names for nml_tracer_data are ’CFC11’ and ’CFC12’ (See Section 13.3)
Required Namelist nml_force_data to specify partial gas pressure (See Section 14.10)
file_data File that contains mole fraction in dry air of CFC11 (and CFC12) [ppt]name CFC11 (and 12)
txyu ’t’

12.3 SF6
The calculation of SF6 also follows the protocols of OMIP (Orr et al., 2016). The formulation of SF6 is nearly the sameas CFCs. Only the coeﬃcients of the empirical fit diﬀer.

12.3.1 Surface boundary condition
See Section 12.2.1, where CFCs should read SF6 .

12.3.2 Usage
Required Model option SF6
Required Add one to numtrc_p
Required Restart file (nmlrs_ptrc) (See Section 11.6)
Required Names for nml_tracer_data is ’SF6’ (See Section 13.3)
Required Namelist nml_force_data to specify partial gas pressure (See Section 14.10)
file_data File that contains atmospheric SF6 [ppt]name SF6
txyu ’t’– 138 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017

Chapter 13

Package Structure and Usage
13.1 Package structure
Packages relevant to tracers are listed as follows.

13.1.1 Tracer equation
tracer_ctl.F90:
tracer.F90:
upc_adv.F90:
quick_adv.F90:
utzq_adv.F90:
som_adv.F90:
mpdata_adv.F90:
tracer_vars.F90:+vvdimp/trcimp.F90:+isopycnal/ipcoef.F90:+isopycnal/ipycmix.F90:

Controller of this package
Main program of this package
Upcurrent advection scheme
QUICK advection scheme (QUICKADVEC)
Combination of UTOPIA and QUICKEST advection scheme (UTZQADVEC)
Second order moment advection scheme (SOMADVEC)MPDATA advection scheme (MPDATAADVEC)
Setting of tracer attributes
Solver of the vertical diﬀusion part using the implicit method (VVDIMP)
Calculation of tensor components of neutral physics parameterization (ISOPYCNAL)
Calculation of tendency due to neutral physics parameterization (ISOPYCNAL)

13.1.2 Vertical mixing coeﬃcientsvmixcoef_ctl.F90:
Controller of the vertical mixing package
vmixcoef.F90:
Main program of the vertical mixing package
vmixcoef_vmbg.F90:
Estimation of background vertical diﬀusion coeﬃcient
vmixcoef_vars.F90:
Declaration of variables+runoﬀ/vmixcoef_rivermouth.F90: Estimation of vertical mixing coeﬃcient around the river mouth

13.1.3 Stratification and convective adjustment
strat_adjust_ctl.F90:
stratification.F90:
cnvajs.F90:
strat_adjust_vars.F90:

Controller of stratification and adjustment package
Main program of calculation of stratification
Main program of convective adjustment
Declaration of variables

13.1.4 Reference state and restoring coeﬃcient
restore_cond_ctl.F90:
restore_cond.F90:
force_data.F90:

Controller of reference state and restoring coeﬃcient
Main program of reference state and restoring coeﬃcientService package that handles external forcing data

13.1.5 Passive tracers
ptrc_ctl.F90:
ptrc.F90:

Controller of passive tracer evolution
Main program of passive tracer evolution that mainly treats surface sources and sinks– 139 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201713.2 Handling the initial state
ptrc_subp.F90:

Sub package of passive tracer evolution that describes internal sources and sinks(called from tracer.F90)+ptrc/cfc.F90, sf6.F90, etc. Sub package of passive tracer evolution that describes specialized processesof particular tracers

13.2 Handling the initial state
How to determine the initial state for temperature and salinity is specified in namelist nml_tracer_run. Parameters arelisted on Table 13.1.
Table13.1 Namelist nml_tracer_run

variable namel_rst_tracer_inunits
logicall_rst_uni_strati

logicalstart_rec_uni_strati
end_rec_uni_strati

integer

description.true. : Read restart files specified by
nmlrs_t and nmlrs_s for the initial condition.
.false.: Start condition depends on thel_rst_uni_strati.false.: Start from 3D-distribution at thestarting time of reference data following
nml_tracer_data.
.true. : Start from uniform stratification
created by the reference data following
nml_tracer_data. Time average is conducted based on start_rec_uni_strati
and end_rec_uni_strati.
uniform stratification is created by the average from start_rec_uni_strati data
record to end_rec_uni_strati record.

usage
Default is the same as l_rst_inof nml_run_ini_state.

if l_rst_tracer_in = .false.

if l_rst_uni_strati = .true.

13.3 Configuration of tracers
The attributes of each tracer such as name, advection scheme, restoring condition, reference data, and restoring coeﬃcients,
are stored in the structural type type_tracer_data defined in tracer_vars.F90. The contents of this structural typeare specified by namelist nml_tracer_data, which should be repeatedly defined as many times as the number of tracersthat should be calculated. Tables 13.2 through 13.8 list the variables.

13.3.1 Name
List of eﬀective names is found in subroutine tracer_vars__set_num_and_name of tracer_data.F90
Table13.2 Namelist nml_tracer_data

variable name
name

units
characterdescription
Name of tracer. Two tracers are necessary:"Potential Temperature" and "Salinity."

usage
Case sensitive. For example, "potential
temperature" is not correct.

13.3.2 Advection scheme
Following can be specified as the name of the advection scheme (adv_scheme%name).
• "upc" : weighted UP-Current advection scheme (always available)• "quick" : QUICK advection scheme (QUICKADVEC)• "utzq" : UTOPIA + ZQUICKEST schemes with ultimate limiter (UTZQADVEC)• "som" : Second-Order Moment advection scheme (SOMADVEC)• "mpdata" : MPDATA advection scheme (MPDATAADVEC)– 140 –TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 2017
Chapter 13

Package Structure and Usage
Table13.3 Namelist nml_tracer_data related to advection scheme

variable nameadv_scheme%nameunitscharacterdescription
Name of the advection
scheme.

adv_scheme%limiter_som_org

logicaladv_scheme%limiter_som_Merryfield03

logicaladv_scheme%lrstin_som

logicaladv_scheme%lrstout_som

logicaladv_scheme%limiter_mpdata_nonoscillatory

logicaladv_scheme%eps_lim_mpdatathat oftracerthat oftracerUse flux limiter for SOM
by Prather (1986)
Use flux limiter for som by
Merryfield and Holloway(2003)The SOM initial state of
moments is read from file
The SOM final state of moments is written to file
Apply flux limiter for MPDATA
Very small value to avoid
zero division
Minimum value for tracer

adv_scheme%min_value_mpdata

usage
Diﬀerent advection
schemes can be set
for individual tracers.
SOMADVEC
SOMADVEC

SOMADVEC
SOMADVECMPDATAADVECMPDATAADVECMPDATAADVEC

13.3.3 Restoring condition
The following are variables related to the restoring condition for a tracer.
Table13.4 Namelist nml_tracer_data related to restoring condition

variable namerestore_conf%l_surf_restore
restore_conf%l_body_restore

units
logicallogicaldescription
restore condition at the surface is applied or not
restore condition in the interior is applied or not

usage
default = .false.
default = .false.

13.3.4 Reference data
a. Three dimensional reference state for restoring

When the field of a tracer is restored to a reference state, the attributes of the reference state should be given by the variables
listed on Table 13.5. This reference state is also used to produce an initial state for that tracer when its restart file is notavailable (See Table 13.1).
Table13.5 Namelist nml_tracer_data related to reference values for body forcing and initial condition

variable nametrcref_conf%file_data

unitscharactertrcref_conf%file_data_grid
trcref_conf%imfrc
trcref_conf%jmfrc
trcref_conf%kmfrc
trcref_conf%interval

character
integer
integer
integer
integer

description
a file name that contains reference
values for body forcing and initialcondition.
a file name of grid
grid size of data in x direction
grid size of data in y direction
grid size of data in z direction
regular time interval of data

– 141 –

usage

needed if linterp = .true.

positive value : unit is sec−1 : monthly−999 : steady forcing
Continued on next page

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201713.3 Configuration of tracers

variable nametrcref_conf%num_data_max
trcref_conf%ifstart

trcref_conf%lrepeat
trcref_conf%linterp
trcref_conf%linterp_v
trcref_conf%ilinear
trcref_conf%luniform
trcref_conf%luniform_v
trcref_conf%ldouble
trcref_conf%iverbose
trcref_conf%ldefined

Table 13.5 – continued from previous page
units
description
integerthe number of record contained in thefileifstart(6) [ymdhms] of the first record of the
input file

logicallogicallogicalinteger
integerlogicallogicalinteger
integerclimatological data is repeatedly usedinterpolate horizontally or notinterpolate vertically or notinterpolation methoddata is horizontally uniform or notdata is vertically uniform or not
input data is double or notstandard output of progress
the input data is defined or not

usage

1999,1,1,0,0,0 when the firstrecord is the average value ofJan 1999 and its data interval
is monthly.
default = .false.
default = .false.
default = .false.
1 : linear, 2 : spline
default = .false.
default = .false.
default = .false.
1 : extensive, 0 : minimum
default = .false.

Format of tracer reference / restoring data is shown in the following.

Format of tracer reference / restoring data (trcref(_surf)_conf%file_data)

integer(4), parameter :: imn = 12, nu = 99
integer(4) :: imfrc, jmftc, kmfrc! data size
character(128) :: file_data, fname_gridreal(4) :: ttlev(imfrc,jmfrc,kmfrc,imn)real(8) :: alonf(imfrc), alatf(jmfrc), dpf(kmfrc)
logical :: linterp, linterp_v! main data
open (unit=nu,file=file_data,access=direct,recl=4*imfrc*jmfrc*kmfrc)
do m = 1, imn
write(unit=nu,rec=m) ttlev(:,:,:,m)end doclose(nu)! longitude/latitude of main data
if (linterp) then ! If input data is horizontally interpolated in the model.
open (unit=nu,file=file_grid)
write(nu) alonf, alatf
if (linterp_v) then ! If input data is vertically interpolated in the model.
write(nu) dpf
end ifclose(nu)end if– 142 –

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

Package Structure and Usageb. Two dimensional reference state for surface restoring

When the surface field of a tracer is intended to be restored to a reference state, the attributes of the surface reference stateshould be given by the variables listed on Table 13.6.
Table13.6 Namelist nml_tracer_data related to reference values for surface restoring forcing.

variable nametrcref_surf_conf%file_data

unitscharactertrcref_surf_conf%file_data_grid
trcref_surf_conf%imfrc
trcref_surf_conf%jmfrc
trcref_surf_conf%interval

character
integer
integer
integertrcref_surf_conf%num_data_max

integertrcref_surf_conf%ifstart

ifstart(6)trcref_surf_conf%lrepeat

logicaltrcref_surf_conf%linterp
trcref_surf_conf%ilinear
trcref_surf_conf%luniform

logicalintegerlogicaltrcref_surf_conf%ldouble
trcref_surf_conf%iverbose
trcref_surf_conf%ldefinedlogicalintegerlogicaldescription
a file name that contains reference values for surface restoring forcing
a file name of grid
grid size of data in x direction
grid size of data in y direction
regular time interval of data

the number of record containedin the file[ymdhms] of the first record of
the input file

climatological data is repeatedly usedinterpolate horizontally or notinterpolation methoddata is horizontally uniform ornotinput data is double or notstandard output of progress
the input data is defined or not

usage

needed if linterp = .true.

positive value : unit is sec−1 : monthly−999 : steady forcing

1999,1,1,0,0,0 when the firstrecord is the average value ofJan 1999 and its data interval
is monthly.
default = .false.
default = .false.
1 : linear, 2 : spline
default = .false.
default = .false.
1 : extensive, 0 : minimum
default = .false.

13.3.5 Restoring coeﬃcienta. Coeﬃcient for three dimensional restoring

When the field of a tracer is restored to a reference state, the attributes of the file that contains restoring coeﬃcients should
be given by the variables listed on Table 13.7. Units of restoring coeﬃcient is sec−1 .
Table13.7 Namelist nml_tracer_data related to restoring coeﬃcient for body forcing.

variable namerstcoef_conf%file_data

unitscharacterrstcoef_conf%file_data_grid
rstcoef_conf%imfrc
rstcoef_conf%jmfrc
rstcoef_conf%kmfrc
rstcoef_conf%interval

character
integer
integer
integer
integer

description
a file name that contains restoring
coeﬃcient for body forcing.
a file name of grid
grid size of data in x direction
grid size of data in y direction
grid size of data in z direction
regular time interval of data

– 143 –

usage
needed if linterp = .true.

positive value : unit is sec
Continued on next page

TECHNICAL REPORTS OF THE METEOROLOGICAL RESEARCH INSTITUTE No.80 201713.3 Configuration of tracers

variable name

Table 13.7 – continued from previous page
units
description

rstcoef_conf%num_data_max

integer

rstcoef_conf%ifstart

ifstart(6)rstcoef_conf%lrepeat

logicalrstcoef_conf%linterp
rstcoef_conf%linterp_v
rstcoef_conf%ilinear
rstcoef_conf%luniform
rstcoef_conf%luniform_v
rstcoef_conf%ldouble
rstcoef_conf%iverbose
rstcoef_conf%ldefinedlogicallogicalintegerlogicallogicallogicallogicallogicalthe number of record contained in
the file[ymdhms] of the first record of the
input file

climatological data is repeatedlyusedinterpolate horizontally or notinterpolate vertically or notinterpolation methoddata is horizontally uniform or notdata is vertically uniform or not
input data is double or notstandard output of progress
the input data is defined or notusage−1 : monthly−999 : steady forcing

1999,1,1,0,0,0 when the firstrecord is the average value ofJan 1999 and its data interval
is monthly.
default = .false.
default = .false.
default = .false.
1 : linear, 2 : spline
default = .true.
default = .false.
default = .true.
1 : extensive 0 : minimum
default = .false.

Note that the default settings for luniform and ldouble are diﬀer from those of the other attributes.
b. Coeﬃcient for surface restoring

When the surface field of a tracer is intended to be restored to a reference state, the attributes of the file that contains surface
restoring coeﬃcients should be given by the variables listed on Table 13.6. Units of the surface restoring coeﬃcient issec−1 .
Table13.8 Namelist nml_tracer_data related to restoring coeﬃcient for surface restoring forcing.

variable namerstcoef_surf_conf%file_data

unitscharacterrstcoef_surf_conf%file_data_grid
rstcoef_surf_conf%imfrc
rstcoef_surf_conf%jmfrc
rstcoef_surf_conf%interval

character
integer
integer
integer

rstcoef_surf_conf%num_data_max
rstcoef_surf_conf%ifstart

ifstart(6)rstcoef_surf_conf%lrepeat

logicalrstcoef_surf_conf%linterp

logicalrstcoef_surf_conf%ilinear
rstcoef_surf_conf%luniform

integerlogicalrstcoef_surf_conf%ldouble
rstcoef_surf_conf%iverbose
rstcoef_surf_conf%ldefinedlogicallogicallogicaldescription
a file name that contains
restoring coeﬃcient for forsurface restoring forcing.
a file name of grid
grid size of data in x direction
grid size of data in y direction
regular time interval of data

the number of record contained in the file[ymdhms] of the first recordof the input file

climatological data is repeatedly usedinterpolate horizontally ornotinterpolation methoddata is horizontally uniform
or notinput data is double or notstandard output of progress
the input data is defined or not– 144 –

usage

needed if linterp = .true.

positive value : unit is sec−1 : monthly−999 : steady forcing

1999,1,1,0,0,0 when the firstrecord is the average value ofJan 1999 and its data interval
is monthly.
default = .false.
default = .false.
1 : linear, 2 : spline
default = .false.
default = .false.
1 : extensive, 0 : minimum
default = .false.

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

Package Structure and Usage

13.3.6 Example
Following is an example of namelist nml_tracer_data for Salinity. Some systems may not allow blank lines or comment
lines in a namelist. In this case, you should delete them.
&nml_tracer_data
name="Salinity",! advection scheme
adv_scheme%name="som",
adv_scheme%limiter_som_org=.false.,
adv_scheme%limiter_som_Merryfield03=.true.,
adv_scheme%lrstin_som=.false.,
adv_scheme%lrstout_som=.true.,! restore_condition
restore_conf%l_surf_restore=.true.
restore_conf%l_body_restore=.false.
! trcreftrcref_conf%file_data=’../data/file_sclim.grd’,trcref_conf%file_data_grid=’dummy.d’,trcref_conf%imfrc=184,
trcref_conf%jmfrc=152,
trcref_conf%kmfrc=51,
trcref_conf%interval=-1,
trcref_conf%ifstart=1947,12,1,0,0,0,
trcref_conf%num_data_max=14,
trcref_conf%lrepeat=.false.,
trcref_conf%linterp=.false.,
trcref_conf%ilinear=1,
trcref_conf%iverbose=1,! rstcoef
rstcoef_conf%ldefined=.false.,! trcref_surf
trcref_surf_conf%file_data=’../data/file_ssurf.grd’,trcref_surf_conf%file_data_grid=’dummy.d’,trcref_surf_conf%imfrc=184,
trcref_surf_conf%jmfrc=152,
trcref_surf_conf%interval=-1,
trcref_surf_conf%ifstart=1947,12,1,0,0,0,
trcref_surf_conf%num_data_max=14,
trcref_surf_conf%lrepeat=.false.,
trcref_surf_conf%linterp=.false.,
trcref_surf_conf%ilinear=1,
trcref_surf_conf%iverbose=1! rstcoef_surf
rstcoef_surf_conf%file_data=’../data/rstcoef_surf_s.grd’,
rstcoef_surf_conf%file_data_grid=’dummy.d’,rstcoef_surf_conf%imfrc=1,
rstcoef_surf_conf%jmfrc=1,
rstcoef_surf_conf%interval=-999,
rstcoef_surf_conf%num_data_max=1,
rstcoef_surf_conf%iverbose=1,
rstcoef_surf_conf%luniform=.true.,
rstcoef_surf_conf%ldouble=.true.,/– 145 –