Description of MIROC6 AGCM

MIROC6 AGCM document writing team∗April 13, 2021

∗ MIROC6 AGCM document writing team (in alphabetical order)

Authors: Taigo Ando, Taro Higuchi, Haruka Hotta, Tomoki Iwakiri, Takuya
Jinno, Kanon Kino, Yuki Takano, Masaki Toda, and Kazuya Yamazaki
Supervisors/Editors: Minoru Chikira, Takanori Kodama, Takuro Michibata, Hiroaki Miura, Tomoko Nitta, Tomoo Ogura, Fuyuki Saito, Miho Sekiguchi, Tatsuo
Suzuki, Kentaroh Suzuki, Hiroaki Tatebe, Masahiro Watanabe, Shingo Watanabe,
and Kei Yoshimura

Remarks

How to cite
MIROC6 AGCM Document Writing Team (2021), Description of MIROC6 AGCM,
CCSR Report No. 65, Division of Climate System Research, Atmosphere and Ocean
Research Institute, The University of Tokyo, https://doi.org/10.15083/0002000180

Publisher Information
Division of Climate System Research,
Atmosphere and Ocean Research Institute,
The University of Tokyo
5-1-5, Kashiwanoha, Kashiwa-shi,
Chiba, 277-8568, Japan
Phone : 04-7136-4371 Fax : 04-7136-4375
Contact: Masahiro Watanabe (hiro@aori.u-tokyo.ac.jp)© 2021 MIROC6 AGCM document writing team1

Acknowledgment
The original ﬁrst version of this document was written in Japanese by Dr. Numaguti (deceased) in 1995 about the former version of the model, CCSR/NIES AGCM5.4, which was
created by himself. Since then, the AGCM has been signiﬁcantly upgraded with continuous community eﬀorts, and it has been periodically renamed by newer version number (e.g.,
CCSR/NIES AGCM5.6, MIROC3, MIROC4, and MIROC5) with almost the same phases of
releases of IPCC’s Assessment Reports (TAR, AR4, and AR5). Even though there are multiple
published papers (e.g., K-1 model developers (2004) for MIROC3 and Watanabe et al. (2010)
for MIROC5) explained about the diﬀerent versions of the model, the description were madeonly for key points without details, and there were no major updates of the original descriptive
document by Dr. Numaguti for a long time.
Here, with the newest release of MIROC series, MIROC6 (Tatebe et al., 2019), it was decided
to make an overhaul of the document for the ﬁrst time since the late 1990’s. This decision was
triggered by some students’ simple question: “Why is this document so obsolete?” Then, under
the TOUGOU project, some of the PIs organized the MIROC6 AGCM document writing teamin the spring of 2020. The team consists of two groups; the authors and the supervisors/editors.
All of the authors are doctor course students of the University of Tokyo, who use or are interested
in using MIROC6, and the supervisors/editors are researchers, who contributed to developmentof MIROC6 and/or its previous versions.
In the beginning, the team converted the original Japanese document written in LaTeX
format into English Markdown format. Ms. Kanon Kino, Ms. Haruka Hotta, Mr. Yuki Takano,and Dr. Fuyuki Saito contributed to make a semi-automated tool for this conversion. Then theteam was divided into several groups, and each group became in charge of each section. These
groups and corresponding sections are as follows:• Mr. Tomoki Iwakiri, Mr. Masaki Toda, and Mr. Kazuya Yamazaki supervised by Dr. Fuyuki
Saito for Dynamics (Ch.2) and Coupler Scheme (Ch.4.1)• Mr. Yuki Takano supervised by Dr. Minoru Chikira for Cumulus Scheme (Ch.3.2)• Mr. Takuya Jinno supervised by Dr. Tomoo Ogura for Shallow Convection Scheme (Ch.3.3)• Ms. Haruka Hotta supervised by Dr. Takuro Michibata and Dr. Kentaroh Suzuki for Large
Scale Condensation (Ch.3.4) and Cloud Microphysics (Ch.3.5)• Mr. Taro Higuchi supervised by Dr. Takanori Kodama and Dr. Miho Sekiguchi for Radiation Scheme (Ch.3.6)• Mr. Taigo Ando supervised by Dr. Minoru Chikira for Trubulence Scheme (Ch.3.7)• Ms. Kanon Kino supervised by Dr. Tatsuo Suzuki for Surface Flux Scheme (Ch.3.8)
It took about six months to draft, and the whole draft was reviewed by all supervisors/editors.
At last, the rest of the sections were amended and improved by all. The team thanks the supportof “Integrated Research Program for Advancing Climate Models (TOUGOU Program)” from the
Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Japan. Finally, theteam would sincerely express our respect and condolences to Dr. Numaguti, the ﬁrst developerof the model and the author of the original version of this document.
Hiroaki Tatebe, Masahiro Watanabe, and Kei Yoshimura
lead member of the teamMarch 31, 2021

2

Revision History
March 31, 2021 First version completed.
April 13, 2021 Add Remarks.

3

Contents
1 Introduction
1.1 Characteristics of MIROC6 AGCM . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Features and Structure of the Model . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Basic Features of the Model. . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Model Execution Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3 Prognostic Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Basic Settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Coordinate System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 Physical Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
8
8
9
10
11
11
12

2 Dynamics
2.1 Basic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Basic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Vertical Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Model Levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Vertical Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Diﬀerences from the σ-Coordinate . . . . . . . . . . . . . . . . . . . . . .
2.3 Horizontal Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Spectral Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Horizontal Diﬀusion Term . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Spectral Representation of Equations . . . . . . . . . . . . . . . . . . . . .
2.4 Time Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Time Integration and Time Filtering with the Leap Frog Method . . . . .
2.4.2 Semi-Implicit Time Integration . . . . . . . . . . . . . . . . . . . . . . . .
2.4.3 Applying the Semi-Implicit Time Integration . . . . . . . . . . . . . . . .
2.4.4 Time Scheme Properties and Requiments for Time Steps . . . . . . . . .
2.4.5 Handling of the Initiation of Time Integration . . . . . . . . . . . . . . . .
2.5 Tracer Advection Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 Introduction of Tracer Advection Scheme . . . . . . . . . . . . . . . . . .
2.5.2 Principle of the Tracer Advection Scheme . . . . . . . . . . . . . . . . . .
2.5.3 The Piecewise Parobolic Method Scheme . . . . . . . . . . . . . . . . . .
2.5.4 Actual Tracer Advection Scheme in MIROC6 . . . . . . . . . . . . . . . .
2.6 Summary of the Dynamical Core . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.1 Conversion of Horizontal Wind to Vorticity and Divergence . . . . . . . .
2.6.2 Calculating a Virtual Temperature . . . . . . . . . . . . . . . . . . . . . .
2.6.3 Calculating the Pressure Gradient Term . . . . . . . . . . . . . . . . . . .
2.6.4 Diagnosis of Vertical Dlow . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.5 Tendency Terms due to Advection . . . . . . . . . . . . . . . . . . . . . .
2.6.6 Transformation of Prognostic Variables to Spectral Space . . . . . . . . .
2.6.7 Transformation of Tendency Terms to Spectral Space . . . . . . . . . . .
2.6.8 Time Integration in Spectral Space . . . . . . . . . . . . . . . . . . . . . .
2.6.9 Transformation of Prognostic Variables to Grid Point Values . . . . . . .
2.6.10 Diﬀusion Correction along Pressure Level . . . . . . . . . . . . . . . . . .
2.6.11 Frictional Heat Associated with Diﬀusion . . . . . . . . . . . . . . . . . .
2.6.12 Horizontal Diﬀusion and Rayleigh Friction . . . . . . . . . . . . . . . . . .
2.6.13 Time Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.14 Correction for Conservation of Mass . . . . . . . . . . . . . . . . . . . . .

14
14
14
15
16
16
17
18
19
19
20
21
24
24
24
25
26
27
28
28
28
29
31
36
36
36
36
37
37
38
38
39
40
41
41
41
42
42

4

2.7

Computational Flow of Dynamical Core . . . . . . . . . . . . . . . . . . . . . . .
2.7.1 Overview of Dynamical Core . . . . . . . . . . . . . . . . . . . . . . . . .

44
44

3 Physics
3.1 Overview of Physical Parameterizations . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Time Integration of Physical Parameterizations . . . . . . . . . . . . . . .
3.1.2 Various Physical Quantities . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Cumulus Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Outline of Cumulus Scheme . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Interaction between Cumulus Ensemble and Large-Scale Environment . .
3.2.3 Cloud Base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.4 Updraft Velocity and Entrainment Rate . . . . . . . . . . . . . . . . . . .
3.2.5 Normalized Mass Flux and Updraft Properties . . . . . . . . . . . . . . .
3.2.6 Spectral Representation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.7 Cloud-Base Mass Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.8 Microphysics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.9 Evaporation, Sublimation and Downdraft . . . . . . . . . . . . . . . . . .
3.2.10 Cloudiness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.11 Cumulus Momentum Transport . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Shallow Convection Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Overview of Shallow Convection . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Basics of Cloud Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Computation in PSHCN . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Large Scale Condensation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Physical Basis for Statistical PDF Scheme . . . . . . . . . . . . . . . . . .
3.4.2 Hybrid Prognostic Cloud scheme . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 PDF Change Through Processes . . . . . . . . . . . . . . . . . . . . . . .
3.4.4 Solving Procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Cloud Microphysics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.1 Overview of Cloud Microphysics . . . . . . . . . . . . . . . . . . . . . . .
3.5.2 Microphysical Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6 Radiation Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6.1 Summary of the Radiation Flux Calculation . . . . . . . . . . . . . . . . .
3.6.2 Wavelength and Sub-Channel . . . . . . . . . . . . . . . . . . . . . . . . .
3.6.3 Calculation of the Planck Function . . . . . . . . . . . . . . . . . . . . . .
3.6.4 Calculation of the Optical Thickness to Gas Absorption . . . . . . . . . .
3.6.5 Calculation of the Optical Thickness to CFC Absorption . . . . . . . . . .
3.6.6 Optical Thickness to Scattering and Scattering Moment . . . . . . . . . .
3.6.7 Total Optical Thickness . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6.8 Expansion of the Plank Function . . . . . . . . . . . . . . . . . . . . . . .
3.6.9 Transmission and Reﬂection Coeﬃcients, and Source Function . . . . . .
3.6.10 T, R, S Matrixes for Maximal/Random Approximation . . . . . . . . . .
3.6.11 Adding of Source Functions for Each Layer . . . . . . . . . . . . . . . . .
3.6.12 Adding in the Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6.13 Calculation of the Temperature Derivative of the Flux . . . . . . . . . . .
3.6.14 Calculation of the Heating Rate . . . . . . . . . . . . . . . . . . . . . . .
3.6.15 Flux of incidence and incident angle . . . . . . . . . . . . . . . . . . . . .
3.6.16 Reading Each Parameter . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6.17 Other Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45
45
45
46
48
48
49
49
50
51
52
53
53
54
54
54
56
56
56
57
61
61
61
62
64
70
70
70
77
77
77
79
79
80
80
84
84
85
88
90
93
94
94
94
96
97

5

3.7

3.8

3.9

Turbulence Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
3.7.1 Surface Layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
3.7.2 Calculation of the Buoyancy Coeﬃcients . . . . . . . . . . . . . . . . . . . 99
3.7.3 Stability Functions for the Level 2 . . . . . . . . . . . . . . . . . . . . . . 100
3.7.4 Master Length Scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
3.7.5 Calculation of Diﬀusion Coeﬃcients . . . . . . . . . . . . . . . . . . . . . 103
3.7.6 Calculation of Turbulent Variables . . . . . . . . . . . . . . . . . . . . . . 106
3.7.7 Time Integration with Implicit Scheme . . . . . . . . . . . . . . . . . . . . 109
Surface Flux Scheme (Sea Surface) . . . . . . . . . . . . . . . . . . . . . . . . . . 111
3.8.1 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
3.8.2 Calculation of Momentum, Heat and Water Vapor Fluxes . . . . . . . . . 115
3.8.3 Radiation Flux Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . 117
3.8.4 Solving Heat Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Gravity Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
3.9.1 Overview of a Gravity Wave Drag Parameterization . . . . . . . . . . . . 120
3.9.2 Relationship between Local FroudeNnumber and Momentum Flux . . . . 120
3.9.3 Momentum Fluxes at the Surface . . . . . . . . . . . . . . . . . . . . . . . 121
3.9.4 Momentum Fluxes in the Upper Levels . . . . . . . . . . . . . . . . . . . 121
3.9.5 The Magnitude of the Time Variation of Horizontal Wind due to Momentum Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
3.9.6 Other Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

4 Miscellaneous
123
4.1 Coupler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
4.1.1 Fluxes to Atmospheric Models . . . . . . . . . . . . . . . . . . . . . . . . 123
4.1.2 Fluxes between Land Surface Model and River Model . . . . . . . . . . . 124
4.1.3 Fluxes to the Ocean Model . . . . . . . . . . . . . . . . . . . . . . . . . . 126
4.2 Deﬁnition of Land-Sea Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . 131
5 References

132

6

1

Introduction

This document describes about the atmospheric general circulation model (AGCM) of MIROC6(the sixth version of MIROC or Model for Interdisciplinary Research on Climate; Tatebe etal., 2019), which has been cooperatively developed at the Center for Climate System Research(CCSR; the precursor of a part of the Atmosphere and Ocean Research Institute), the University of Tokyo, the Japan Agency for Marine-Earth Science and Technology (JAMSTEC), andthe National Institute for Environmental Studies (NIES). All the descriptions are basically corresponded with the source codes of MIROC6, speciﬁcally the version used for CMIP6 DECK
experiment, whereas the original document written in Japanese by Dr. Numaguti in 1995 wascorresponded with those of CCSR/NIES AGCM5.4.

1.1

Characteristics of MIROC6 AGCM

MIROC6 AGCM are summarized below.
• System of equations: Hydrostatic primitive equations• Area: Global 3D• Prognostic variables: Horizontal wind speed, temperature, surface pressure, speciﬁc
humidity, cloud water• Horizontal discretization: Spectral transformation (Bourke, 1988) method• Vertical discretization: Hybrid σ − p coordinate, based on Arakawa and Konor (1996)• Resolution for default: T85 (150 km), 81 levels up to 0.004 hPa• Time integration: Essentially the leap frog scheme, with a time ﬁlter (Williams, 2009)• Cumulus: An entrainment plume model with multiple cloud types (Chikira and Sugiyama,
2010)• Shallow convection: A mass-ﬂux-based single-plume model based on Park and Bretherton (2009)• Large scale condensation & Cloud microphysics: A prognostic large scale condensation scheme (Watanabe et al., 2009) and the implementation of a bulk mirco-physicalscheme (Wilson and Ballard, 1999)• Radiation: k-distribution scheme (Sekiguchi and Nakajima, 2008) with a hexagonal solid
column as ice particle habit and extended mode radius of cloud particles• Turbulence: The Mellor-Yamada-Nakanishi-Niino scheme (Nakanishi 2001; Nakanishiand Niino 2004)’s level 2.5 closure scheme• Surface ﬂux: Bulk coeﬃcients (Louis, 1979; Louis et al., 1982) with convection eﬀects at
sea surface (Miller et al., 1992)• Gravity wave drag: An orographic gravity wave parameterization (McFarlane, 1987)
with a non-orographic gravity wave parameterization (Hines, 1997; Watanabe et al., 2011)

7

1.2
1.2.1

Features and Structure of the Model
Basic Features of the Model.

The MIROC6 AGCM is a numerical model for describing the global three-dimensional atmosphere based on physical laws and calculating the time evolution of the system as an initial value
problem or a boundary value problem.
The data to be inputted are as follows.
• Initial data for each prognostic variable (horizontal wind speed, temperature, surface pressure, speciﬁc humidity, cloud liquid water content, etc.)• Boundary condition data (surface elevation, surface condition, sea surface temperature,etc.)• Various parameters of the model (atmospheric components, physical process parameters,etc.)
On the other hand, the output is the following.
• Data for each prognostic parameter and diagnostic parameter, for each time or time average• Initial data to be used for continuous execution (restart data)• Progress and various diagnostic messages
The prognostic variable is the data obtained as a time series by integrating the diﬀerential
equation of time evolution, and the diagnostic variable is the quantity calculated from theprognostic variable, the boundary conditions and the parameters by some method that does notinclude time integration.
More speciﬁcally, the model basically solves the following equations (prognostic equations).
∂u= (Fx )D + (Fx )P .
∂t∂v
= (Fy )D + (Fy )P .
∂t∂T
= (Q)D + (Q)P .
∂t∂pS= (M )D + (M )P .
∂t∂q
= (S)D + (S)P .
∂t∂Tg= (Qg )D + (Qg )P .
∂t

(1.1)
(1.2)
(1.3)
(1.4)
(1.5)
(1.6)

Here, u, v, T, pS , q, Tg are two-dimensional and three-dimensional prognostic variables suchas eastward wind, northward wind, temperature, surface pressure, speciﬁc humidity, and surface
state amount, respectively, and the right-hand side is a term that causes time variation of eachprognostic variable. The terms Fx , Fy , Q, S, Qg are calculated based on the prognostic variablesu, v, T, pS , q, Tg , are divided into two main categories: the terms u and v, such as advection due
to the motion of the atmosphere (the terms with index D in the above equation), and the terms
such as cloud and radiation (the terms with index P in the above equation). There are two8main types of terms. The former is called the dynamical process, and the latter is called thephysical process.
The advection term is the main part of the time-varying term in dynamical processes, andthe accurate estimation of the spatial derivative is important in its calculation. The MIROC6
AGCM utilizes the spherical harmonic expansion to calculate the horizontal diﬀerential term.
On the other hand, it is important for physical processes to be represented in a simple modelwith parameters (parameterization), such as energy conversions due to the phase change ofwater, radiative absorption and emission, the eﬀects of small-scale atmospheric motions, andthe eﬀects of various processes on the ground surface.
The time integration of the prognostic equation is done by approximating the left-hand sideof (1) etc. by the diﬀerence. For example,∂q
q t+∆t − q t
→
∂t
∆t

(1.7)

q t+∆t = q t + ∆t [(S)D + (S)P ](1.8)

By making,where S is a function of the prognostic variables u, v, T, pS , q. Although S is a function of
the prognostic variables u, v, T, pS , q, and so on, there are various time diﬀerence schemes thatcan be used in this calculation depending on the time of day the prognostic variables are used
to evaluate S. The MIROC6 AGCM uses the Euler method, which uses the value of the t as itis, the leap frog method, which uses the value of the t + ∆t/2, and the implicit method, whichuses the (approximate) value of the t + ∆t.
In the MIROC6 AGCM, the time integration of the prognostic variables is done separately
for the dynamical and physical processes. The dynamical processes basically use a leap frog,q̃ t+∆t = q t−∆t + 2∆t (S)tD(1.9)

However, some terms are treated as implicit. In the physical process, based on the results ofintegrating the dynamical terms, the Euler and implicit methods are used together,q t+∆t = q̃ t+∆t + 2∆t (S)P

(1.10)in (8). Note that ∆t in (8) is replaced by 2∆t.
1.2.2

Model Execution Flow

The ﬂow of the model execution is brieﬂy shown below. The entries in the index are the namesof the corresponding subroutine.
1. set the parameters of an experiment, coordinates, etc.
SUBROUTINE:[PCONST,ASETCO,SETPAR,SETTSTRT,SETTEND]
2. read the initial values of the prognostic variables SUBROUTINE:[RDSTRT]
3. start the time step SUBROUTINE:[TIMSTP]
94. perform time integration by mechanical processes SUBROUTINE:[DYNMCS]
5. perform time integration by physical processes SUBROUTINE:[PHYSCS]
6. advance the time MODULE:[TFILT]
7. Output the data if necessary MODULE:[HISTOU]
8. Output the restart data if necessary SUBROUTINE:[WRRSTR]
9. Return to 3
1.2.3

Prognostic Variables

The prognostic variables are as follows. The values in parentheses are the coordinate system,and λ, φ, σ, z indicate the longitude, latitude, dimensionless pressure, σ, and vertical depth,
respectively. The values in the square brackets are in units of the index.
Element

SymbolUnit

Eastward wind speed
Northward wind speed
Atmospheric temperature
Surface pressure
Speciﬁc humidity
Cloud water speciﬁc humidity
Cloud ice speciﬁc humidity
Total water PDF varianceTotal water PDF skewness
Variance of liquid potential temperature
Covariance of liquid potential temperature and total water
Variance of total water
Tracersu (λ, φ, σ)v (λ, φ, σ)T (λ, φ, σ)pS (λ, φ)q (λ, φ, σ)l (λ, φ, σ)qi (λ, φ, σ)V (λ, φ, σ)S (λ, φ, σ)
T SQ (λ, φ, σ)
COV (λ, φ, σ)
QSQ (λ, φ, σ)
-

[m/s]
[m/s]
[K]
[hPa]
[kg/kg][kg/kg][kg/kg]
[ND]
[ND]
[K2 ]
[K]
[ND]
-

Of these quantities, the quantities for turbulence process, T SQ, COV, QSQ, store only one
step at a time, while the quantities for the atmosphere, u, v, T, pS , q, l, qi , V, S, need to store two
steps at a time. This is due to the fact that the leap frog method is used in the time integrationof the dynamic process of the quantities related to the atmosphere.
The quantities of the atmosphere, u, v, T, pS , q, l, are variables managed by the main routine,
Administration of the Atmosphere'[AGCM5\a]. On the other hand, the quantities relating
to the earth’s surface and ground, qi , V, S, T SQ, COV, QSQ, do not appear in the main routine,
but are managed by the subroutine MODULE:[PHYSCS] of the physical process.
Tracers include mass concentrations of aerosol species,

10

1.3

Basic Settings

Here we present the basic setup of the model.
1.3.1

Coordinate System

The coordinate system of the atmospheric model consists of longitude λ, latitude φ, and normalized pressure η (deﬁnitions are given below), each of which is treated as orthogonal. However, z
is used for the vertical coordinate in the ground, which is treated in a land physics component.
Longitude is discretized at equal intervals (SUBROUTINE: [SETLO] in asetc.F).
λi = 2πi−1, i = 1, . . . I.
I

(1.11)

Latitude grids φj are derived from the Gauss-Legendre integral formula (SUBROUTINE: [SETLA]in asetc.F). This is the zero point of the Legendre polynomial of order J with µ = sin φ as theargument (SUBROUTINE: [GAUSS] in uspst.F). If J is large, we can approximate
1 j − 1/2−
, j = 1, . . . J.
(1.12)
φj = π2
J
Usually, the grid spacing of longitude and latitude is taken to be approximately equal to J = I/2,
based on the triangular truncation of the spectral method.
Air pressure p is deﬁned at half levels (pk+1/2 , k = 1, 2, . . . K) using the following formulausing constants Ak+1/2 , Bk+1/2 :pk+1/2 = Ak+1/2 + Bk+1/2 ps ,(1.13)

where A1/2 = AK+1/2 = 0, B1/2 = 1, BK+1/2 = 0 and thus p1/2 = ps , pK+1/2 = 0. Therefore,
the normalized pressure σ ≡ p/ps can be written as below:
σk+1/2 =Ak+1/2+ Bk+1/2 .
ps(1.14)

Furthermore, a hybrid normalized pressure η is deﬁned as below:ηk+1/2 =Ak+1/2+ Bk+1/2 ,p0p0 ≡ 1000 hPa.

(1.15)

Since Ak+1/2 , Bk+1/2 , p0 are constants, ηk+1/2 is also a constant and we use it as the verticalcoordinate of the atmopheric model. However, as described in Chapter 2, basic equations aredescretized in such a way that ηk+1/2 does not explicitly appear and σk+1/2 is used instead tocommonize source codes with the σ-coordinate system used in MIROC 5.
Air pressure pk at full levels (pk , k = 1, 2, . . . K) is interpolated from half-level pressure as
below:(
κ+1 !)1/κ
pκ+11
k−1/2 − pk+1/2
pk =.
(1.16)
1 + κ pk−1/2 − pk+1/2
Full-level pressure in a 80-level conﬁguration is shown in Fig. 1. While lower layers follow theterrain, upper layers are isobaric, and the two are smoothly connected.
All prognostic variables are deﬁned either on a grid of (λi , φj , ηk ) or (λi , φj , zl ). (The underground level, zl , is described in the section on physical processes.)In the time direction, the forecast equations are discretized at evenly spaced ∆t and time
integration is performed. However, ∆t may change in cases where the stability of the time
integration is insuﬃcient.
11

k=80
10 2

Pressure (hPa)

10 1
100
101
102
103

k=70
k=60
k=50
k=40
k=30
k=20
k=10

Figure 1: Default arangement of vertical levels for 80-level simulations.
1.3.2

Physical Constants

The basic physical constants are shown below (SUBROUTINE [PCONST] in apcon.F).
Element

Symbol

Unit

Value

Earth radius
Gravitational acceleration
Atmospheric speciﬁc heat at constant pressure
Atmospheric gas constantLatent heat of water evaporation
Water vapor speciﬁc heat at constant pressure
Gas constant of water
Density of liquid water
Saturated vapor pressure at 0 ◦ C
Stefan-Bolzman constantKárman constantLatent heat of ice melting
Freezing point of water
Constant pressure speciﬁc heat of water
Freezing point of seawaterSpeciﬁc heat ratio of ice at constant pressure
Water vapor molecular weight ratio
Coeﬃcient of virtual temperature

a
g
Cp
R
L
Cv
Rv
dH2 O
e∗ (273K)
σSB
k
LM
TM
Cw
TI
CI = Cw − LM /TM
ϵ = R/Rv
ϵv = ϵ−1 − 1m
m s−2
J kg−1 K−1
J kg−1 K−1
J kg−1J kg−1 K−1
J kg−1 K−1
kg m−3Pa
W m−2 K−4

6.37 ×1069.8
1004.6287.04
2.5 ×106
1810
461
1000
611
5.67 ×10−8
0.4
3.4 ×105
273.154200
271.35
2397
0.622
0.606

12

J kg−1K
J kg−1K

Element

SymbolRatio of speciﬁc heat to gas constantκ = R/Cp

13

Unit

Value
0.286

2

Dynamics

2.1

Basic Equations

2.1.1

Basic Equations

The basic equations are a system of primitive equations at the spherical (λ, φ) and η coordinates,
given as follows (Arakawa and Konor 1996).
1. Continuity equation∂m∂(mη̇)+ ∇η · (mvH ) +
=0
∂t
∂η

(2.1)

∂Φ
RTv
=−
m
∂η
p

(2.2)

∂ζ
1 ∂Av
1
∂
=
−
(Au cos φ) − D(ζ)∂ta cos φ ∂λa cos φ ∂φ

(2.3)

∂D
1 ∂Au
1
∂
=
+
(Av cos φ) − ∇2η (Φ + RT̄ π + E) − D(D)∂ta cos φ ∂λa cos φ ∂φ

(2.4)

2. Hydrostatic equation

3. Equation of motion

4. Thermodynamic equation∂T∂t

1 ∂uT ′
1 ∂−(vT ′ cos φ) + T ′ D
a cos φ ∂λ
a ∂φ
 


Qdif f
∂π
∂T
κT
mη̇
Q
− η̇
+
B
+ vH · ∇ η π +++− D(T )
∂η
σ
∂t
ps
Cp
Cp

= −

(2.5)

5. Tracers
For any tracer whose conservative quantity (e.g. mixing ratio) is denoted as q,
∂q
∂t

∂
1 ∂uq
1
−
(vq cos φ) + qD
a cos φ ∂λa cos φ ∂φ∂q
− η̇
+ Sq − D(q)
∂η

= −

14

(2.6)

Here,

m ≡

∂p
∂η


,

(2.7)

ps

θ ≡ T (p/p0 )−κ ,(2.8)

κ ≡ R/Cp ,(2.9)

Φ ≡ gz,(2.10)π ≡ ln pS ,
dη
η̇ ≡
,
dt
Tv ≡ T (1 + ϵv q),(2.11)(2.12)(2.13)′

T≡ T̄ + T ,(2.14)

T̄

≡ 300 K,1∂
1 ∂v−
(u cos φ),≡a cos φ ∂λ a cos φ ∂φ
1 ∂u
1
∂
≡
+
(v cos φ),a cos φ ∂λ a cos φ ∂φ∂u
RT ′ ∂π≡ (ζ + f )v − η̇−+ Fx ,∂η a cos φ ∂λ∂v RT ′ ∂π
≡ −(ζ + f )u − η̇−+ Fy ,∂ηa ∂φu2 + v 2
,
≡
2  
 
u
∂
v ∂
≡
+
,
a cos φ ∂λ σ a ∂φ σ


1
∂2
1
∂
∂
≡
+
cos φ
.
a2 cos2 φ ∂λ2 a2 cos φ ∂φ
∂φ(2.15)

ζDAu
Av
E
vH · ∇∇2η

(2.16)(2.17)(2.18)(2.19)(2.20)(2.21)(2.22)

f is the Coriolis parameter. D(ζ), D(D), D(T ), D(q) are horizontal diﬀusion terms, Fλ , Fφ
are forces due to small-scale kinetic processes (treated as ‘physical processes’), Q are forces due
to radiation, condensation, small-scale kinetic processes, etc. Heating and temperature change
due to ‘physical processes’, and Sq is a water vapor source term due to ‘physical processes’ suchas condensation and small-scale motion. Qdiﬀ is the heat of friction and ∂v
Qdiﬀ = −v ·.
(2.23)∂t diﬀ( ∂v∂t )diﬀ is a time-varying term of u, v due to horizontal and vertical diﬀusion.
2.1.2

Boundary Conditions

Upper and lower boundary conditions for the vertical velocity is:η̇ = 0atη = 0, 1.

(2.24)The prognostic equation for ps and the diagnostic equation for the vertical velocity can bederived by integrating the continuity equation and applying these boundary conditions.

15

2.2

Vertical Discretization

Following Arakawa and Konor (1996) except for using the Lorentz grid, the basic equations are
discretized vertically by diﬀerences. This scheme has the following characteristics.
• Conservation of the total integrated mass.
• Conservation of the total integrated energy.
• Conservation of the globally integrated angular momentum.
• Conservation of the total mass-integrated potential temperature.
• The hydrostatic pressure equation is localized (the altitude of the lower level is independentof the temperature of the upper level).
• For a given temperature distribution, constant in the horizontal direction, the hydrostaticpressure equation becomes precisely accurate and the barometric gradient force becomeszero.
• Isothermal atmosphere stays permanently isothermal.
2.2.1

Model Levels

Model level increases in altitude with the vertical level number k. k = 1/2 corresponds withthe model bottom (η = 1), while k = K + 1/2 corresponds to the model top (η = 0). Variables
ζ, D, T, q are deﬁned at integer levels (k = 1, 2, . . . K), while the vertical velocity η̇ is deﬁned at
half-integer levels (k = 1/2, 3/2, . . . K + 1/2). Using constants Ak+1/2 and Bk+1/2 and variablesurface pressure ps , air pressure at half levels are deﬁned as below:
pk+1/2 = Ak+1/2 + Bk+1/2 ps .

(2.25)

Thus, the normalized pressure σ ≡ p/ps can be written as below:
σk+1/2 =Ak+1/2+ Bk+1/2 .
ps(2.26)

Using a reference pressure p0 , the hybrid-normalized pressure η is deﬁned as below:ηk+1/2 =Ak+1/2+ Bk+1/2 ,p0(2.27)

which is a constant at all levels and is used as the vertical coordinate by default in MIROC 6.0.
Pressure at full levels are interpolated from half-level pressure by the following formula:
(
pk =

1
1+κ

κ+1
pκ+1
k−1/2 − pk+1/2!)1/κ

pk−1/2 − pk+1/2.

(2.28)For later use, let us deﬁne the following:∆σk ≡ σk−1/2 − σk+1/2 ,∆Bk ≡ Bk−1/2 − Bk+1/2 .

16(2.29)

2.2.2

Vertical Discretization

Basic equations vertically discretized at the η hybrid coordinates are shown below.
1. Continuity equation and diagnosis of the vertical velocityX∂π{Dk ∆σk + (vk · ∇π)∆Bk }=−∂t
K(2.30)

k=1

In MIRCO 6.0, the discretization is conducted in a manner similar to the σ coordinate,
which can be optionally selected and was the default in previous versions, to commonize source
codes. Thus, the vertical velocity is represented as σ̇ = mη̇/ps . Furthermore, vertical advectionη̇(∂/∂η) is replaced with an equivalent form mη̇/ps (∂/∂σ).
(σ̇ =)

K(mη̇)k−1/2∂π X= −Bk−1/2−{Dl ∆σl + (vl · ∇π)∆Bl }ps∂t(2.31)

l=k(mη̇)1/2(mη̇)k+1/2
=
=0
ps
ps

(2.32)

2. Hydrostatic equation
Φ1 = Φs + Cp (σ1−κ − 1)Tv,1= Φs + Cp α1 Tv,1Φk − Φk−1





 
pk−1/2 κ
pk−1/2 κ= Cp− 1 Tv,k + Cp 1 −Tv,k−1pk
pk−1= Cp αk Tv,k + Cp βk−1 Tv,k−1(2.33)

Here,



αk
βk


pk−1/2 κ
≡
− 1,
pk


pk+1/2 κ≡ 1−.
pk(2.34)(2.35)(2.36)

3. Equations of motion
∂ζk
∂t
∂D
∂t

=
=

(Au )k =(Av )k =κ̂k =1 ∂(Av )k1∂−(Au cos φ)k − D(ζk )(2.37)
a cos φ ∂λa cos φ ∂φ∂
1 ∂(Au )k1+(Av cos φ)k − ∇2η (Φk + RT̄ π + (K E)k ) − D(Dk(2.38))
a cos φ ∂λa cos φ ∂φ(mη̇)k−1/2 uk−1 − uk(mη̇)k+1/2 uk − uk+1
(ζk + f )vk −+ps∆σk−1 + ∆σkps∆σk + ∆σk+1
1 ∂π(Cp Tv,k κ̂ − RT̄ ) + Fx
(2.39)−
a cos φ ∂λ(mη̇)k−1/2 vk−1 − vk(mη̇)k+1/2 vk − vk+1−(ζk + f )uk −+ps∆σk−1 + ∆σkps∆σk + ∆σk+1
1 ∂π−(Cp Tv,k κ̂ − RT̄ ) + Fy
(2.40)a ∂φ
Bk−1/2 αk + Bk+1/2 βk(2.41)∆σk17

4. Thermodynamic equation(Qdif f )k∂Tk
1 ∂uk Tk′
1
∂
Qk
=−
−
(vk Tk′ cos φ) + Hk ++− D(Tk )∂ta cos φ ∂λa cos φ ∂φCp
Cp

(2.42)

Here,
"

Hk

(mη̇)k−1/2 T̂k−1/2 − Tk(mη̇)k+1/2 Tk − T̂k+1/2≡ Tk′ Dk −
+
ps
∆σk
ps
∆σk
( "
#
K
X
+ αk Bk−1/2 vk · ∇π −(Dl ∆σl + (vl · ∇π)∆Bl )l=k"+ βk Bk+1/2 vk · ∇π −K
X#)
(Dl ∆σl + (vl · ∇π)∆Bl )

l=k+1

"

#

1
Tv,k
∆σk
#

(mη̇)k−1/2 T̂k−1/2 − Tk(mη̇)k+1/2 Tk − T̂k+1/2
+
ps
∆σl
ps
∆σl

= Tk′ Dk −+κ̂k (vk · ∇π)Tv,k
K
X
Tv,k
−αk
(Dl ∆σl + (vl · ∇π)∆Bl )
∆σk

−βk

l=k
K
X

(Dl ∆σl + (vl · ∇π)∆Bl )

l=k+1

Tv,k
,
∆σk

T̂k−1/2 = ak Tk + bk−1 Tk−1 ,


 −1
pk κ
ak = αk 1 −
,
pk−1


−1
pk κ
bk = βk−1.
pk+1(2.43)(2.44)(2.45)(2.46)

5. Tracers
∂qk
∂t
Rk
2.2.3

1 ∂uk qk1∂−
(vk qk cos φ) + Rk + Sq,k − D(qk )a cos φ ∂λa cos φ ∂φ(mη̇)k+1/2 qk − qk+1
1 (mη̇)k−1/2 qk−1 − qk
= qk D k −
+
2
ps
∆σk
ps
∆σk
= −

(2.47)(2.48)

Diﬀerences from the σ-Coordinate

In MIROC 6.0, the discretization is conducted in a similar form to the σ coordinate. Thus,diﬀerences of discretized equations between the η and σ coordinates are relatively small, which
are listed below:• In the σ coordinate, Ak+1/2 is equal to zero at all levels.
• While ∆Bk and ∆σk are diﬀerent in the η coordinates, those are equivalent to each otherin the σ coordinate.

18

2.3

Horizontal Discretization

The horizontal discretization is based on the spectral transformation method (Bourke, 1988).
The diﬀerential terms for longitude and latitude are evaluated by the orthogonal function expansion, while the non-linear terms are calculated on the grid.
2.3.1

Spectral Expansion

As an expansion function, the spherical harmonic functions Ynm (λ, µ), which are eigenfunction of
Laplacian on a sphere, are used. However, µ ≡ sin φ is used. Ynm satisﬁes the following equation,∇2σ Ynm (λ, µ) = −n(n + 1) m
Yn (λ, µ)a2(2.49)

Using the Associated Legendre function Pnm it is written as follows.
Ynm (λ, µ) = Pnm (µ)eimλ(2.50)

However, it is n ≥ |m|.
The expansion by spherical harmonic functions is ,Ynm ij ≡ Ynm (λi , µj )(2.51)

When I write ,NX

Xij ≡ X(λi , µj ) = R⌉NX

Xnm Ynm ij ,(2.52)

m=−N n=|m|The inverse of that is ,

Xnm =

1
4π

Z 1

Z π
dλX(λ, µ)Ynm∗ (λ, µ)

dµ
−1

0

1 XX
=
Xij Ynm∗ ij wjIIJ(2.53)

i=1 j=1

The formula is expressed as To evaluate by replacing the integral with a sum, we use theGauss trapezoidal formula for the λ integral and the Gauss-Legendre integral formula for the µ
integral. µj is the Gauss latitude and wj is the Gaussian weights. Also, λi is a grid of evenly
spaced Gaussian weights.
Using spectral harmonics transformation, the grid point values of the terms containing thederivatives can be calculated as follows.


∂X
∂λ


= Re
ij

N
X

N
X

m=−N n=|m|19

imXnm Ynm ij(2.54)



N
X
∂X
cos φ
= Re
∂φ ij

N
X

Xnm (1 − µ2 )m=−N n=|m|∂ m
Y∂µ n ij(2.55)

Furthermore, the grid point values of u, v can be obtained from the spectral components ofζ and D as follows

N 
N
X
X
a
1
∂
ima
Re
ζnm (1 − µ2 ) Ynm ij −Dnm Ynm ijcos φn(n + 1)∂µn(n + 1)(2.56)


N
N 
X
X
1
ima
a
m m
m
2 ∂
m
vij =
Re
−
ζ Y
−
D (1 − µ ) Yn ij
cos φn(n + 1) n n ij n(n + 1) n∂µ

(2.57)

uij =

m=−N n=|m|n̸=0m=−N n=|m|n̸=0

The derivative appearing in the advection term of the equation is calculated as


m

Z 1
Z π
1
1 ∂A m∗
=
dµ
dλ
Y
4π
a
cos
φ ∂λ n
−1
0
n
Z 1
Z π
1
1
=
dµ
dλ imA cos φ
Y m∗
2) n
4π −1
a(1
−
µ
0

1 ∂A
a cos φ ∂λwj
1 XX
imAij cos φj Ynm∗ ijIa(1 − µ2j )I=

J(2.58)

i=1 j=1



m
Z 1
Z π
1
∂
1
∂
(A cos φ)=dµdλ (A cos φ)Ynm∗a cos φ ∂φ
4πa −1
∂µ
0
n
Z 1
Z π
1
∂
=−
dµ
dλA cos φ Ynm∗
4πa −1
∂µ
0
wj
1 XX
∂
Aij cos φj (1 − µ2j ) Ynm∗ ijI∂µa(1 − µ2j )
I

=−

J

(2.59)

i=1 j=1

Furthermore,mn(n + 1) m
∇2σ X n = −Xna2(2.60)

to be used for the evaluation of the ∇2 section.
2.3.2

Horizontal Diﬀusion Term

The horizontal diﬀusion term is implemented in the form ∇ND as follows.
"D(ζ) = KM H (−1)ND /2 ∇ND −
20

2
a2

ND /2 #ζ,(2.61)

"

D(D) = KM H (−1)

ND /2

∇

ND

−

2
a2

ND /2 #D,(2.62)

D(T ) = (−1)ND /2 KHH ∇ND T,(2.63)

D(q) = (−1)ND /2 KEH ∇ND q.

(2.64)

This horizontal diﬀusion term damps high frequency component occuring aliasing for computational stability. In order to represent selective horizontal diﬀusion on small scales, 4 ∼ 16 isused as ND . Here, the extra term for vorticity and divergence diﬀusion indicates that the termof rigid body rotation in n = 1 does not decay.
2.3.3

Spectral Representation of Equations

1. A series of equations

m
X
∂πm
=−
(Dnm )k ∆σk
∂t
K

+

1
I

k=1
I
J
XX

Zij Ynm∗ ij wj ,(2.65)

i=1 j=1

where,Z≡−K
X

vk · ∇π.

(2.66)

k=1

2. Equation of motion

wj
∂ζnm
1 XX
=
im(Av )ij cos φj Ynm∗ ij∂tIa(1 − µ2j )I

J

i=1 j=1

wj
1 XX∂(Au )ij cos φj (1 − µ2j ) Ynm∗ ijI∂µa(1 − µ2j )I+

J

i=1 j=1m−(DM )mn ζn ,

21(2.67)

wj
1 XX∂ D̃nm=im(Au )ij cos φj Ynm∗ ij∂tIa(1 − µ2j )I

J

i=1 j=1

wj∂
1 XX(Av )ij cos φj (1 − µ2j ) Ynm∗ ijI∂µa(1 − µ2j )I−

J

i=1 j=1n(n + 1) 1 X X
Eij Ynm∗ ij wj
a2
I
I

−

J

i=1 j=1+n(n + 1) mm
(Φn + Cp κ̂k T̄k πnm ) − (DM )mn Dn ,
a2(2.68)where"(DM )mn = KM Hn(n + 1)a2

ND /2−2a2

ND /2 #.

(2.69)

3. Thermodynamic equation

wj
1 XX
∂Tnm
=−
imuij Tij′ cos φj Ynm∗ ij∂tIa(1 − µ2j )I

J

i=1 j=1

+

1
I

I X
J
X

vij Tij′ cos φj (1 − µ2j )i=1 j=1

wj
∂ m∗
Yn ij
∂µ
a(1 − µ2j )IJ 
Qij + Qdif f
1 XX+
Ynm∗ ij wj
Hij +
I
Cp
i=1 j=1m−(D̃H )mn Tn ,(2.70)

where,(DH )mn = KHHn(n + 1)a2

ND /2.

(2.71)

4. Water vapor formula

wj
∂qnm
1 XX
=−
imuij qij cos φj Ynm∗ ij∂tIa(1 − µ2j )I

J

i=1 j=1

+

1
I

I X
J
X

vij qij cos φj (1 − µ2j )i=1 j=1

wj
∂ m∗
Yn ij
∂µ
a(1 − µ2j )


1 XX
+
R̂ij + Sq,ij Ynm∗ ij wj
IIJ

i=1 j=1m+(DH )mn qn

22(2.72)

where(DE )mn = KEHn(n + 1)a2

23

ND /2.

(2.73)

2.4

Time Integration

The time discretization is essentially the leap frog scheme. However, backward or forward
diﬀerences are used for diﬀusion terms and physical process terms. A time ﬁlter (Williams,
2009), which is a modiﬁed version of the Asselin time ﬁlter (Asselin 1972), is used to suppress
computational modes. A semi-implicit method is applied to the gravitational wave term to makethe ∆t larger (Bourke, 1988).
2.4.1

Time Integration and Time Filtering with the Leap Frog Method

We use leap frog as the time integration scheme for advection terms and other dynamic terms.
A backward diﬀerence of 2∆t is used for the horizontal diﬀusion term. The p-surface correctionof the diﬀusion term and the frictional heat due to horizontal diﬀusion are treated by forward
diﬀerences of 2∆t. The physical process terms (Fλ , Fφ , Q, Sq ) use the forward diﬀerence of2∆t (except for the vertical diﬀusion term, which uses the forward diﬀerence of Fλ , Fφ , Q, Sq ).
However, the calculation of the prognositc veriables of vertical diﬀusion is treated as a backward
diﬀerence. Please refer to the chapter on physical processes for details.)
Expressing each prognostic variable as X,
(2.74)
X̂ t+∆t = X̄ t−∆t + 2∆tẊadv X t + 2∆tẊdif X̂ t+∆t ,where Ẋadv is the advection term etc., and Ẋdif is the horizontal diﬀusion term.
X̂ t+∆t is then corrected for diﬀusion (Ẋdis for p-surface correction and the heat of friction)and physical processes (Ẋphy ), yielding X t+∆t .
(2.75)X t+∆t = X̂ t+∆t + 2∆tẊdis X̂ t+∆t + 2∆tẊphy X̂ t+∆t
To damp numerical modes, a time ﬁlter (Williams, 2009) is applied to leap-frog method atevery steps. The time ﬁlter is given below, where terms with over bars are ﬁltered.
¯ t = X̄ t + να[X̄¯ t−∆t − X̄ t + X t+∆t ],X̄

(2.76)¯ t−∆t − 2X̄ t + X t+∆t ],X̄ t+∆t = X t+∆t + ν(1 − α)[X̄(2.77)

where ν = 0.05 and α = 0.5.
2.4.2

Semi-Implicit Time Integration

Basically, the leap frog is used for the dynamic processes, but the trapezoidal implicit scheme
is used for some terms. For a vector quantity q, let us write the value at t as q, the value att + ∆t as q+ , and the value at t − ∆t as q− . Then, in the trapezoidal implicit scheme, the timechange term is evaluated for (q+ + q− )/2, instead of q used in the simple leap forg method. We
now divide q into two time varying terms, one (A) for the leap forg method and the other (B)for the trapezoidal implicit method. We assume that (A) is nonlinear to q, while (B) is linear.
In other words,
q+ = q− + 2∆tA(q) + 2∆tB(q+ + q− )/2,

(2.78)

where (B) is a square matrix. Deﬁning ∆q ≡ q+ − q, we get(I − ∆tB)∆q = 2∆t (A(q) + Bq) .
This can be easily solved by matrix operations.
24(2.79)

2.4.3

Applying the Semi-Implicit Time Integration

Here, we apply the semi-implicit method and treat terms associated with linear gravity waves
as implicit, which allows us to increase the time step ∆t.
We divide the basic equation into a linear gravity wave term (T = T̄k ) with a static ﬁeld
as the basic ﬁeld and other terms (with the indices N G). Using a vector representation for thevertical direction (D = {Dk } and T = {Tk }),
 
∂π
∂π
=
− C · D,(2.80)∂t∂t N G∂D∂D=− ∇2η (ΦS + W T + Gπ) − DM D,(2.81)∂t∂t N G∂T∂T=− hD − DH T.
(2.82)∂t∂t N G
Here, the non-gravitational wave term is


∂π
∂t

N G
= −

K
X

vk · ∇π∆Bk ,

k=1



G
(mη̇)N
k−1/2

∂π
∂t

N G

(2.83)KX= −Bk−1/2−
vl · ∇π∆Bl ,
ps
l=k


∂D N G
1 ∂(Au )k1∂=+(Av cos φ)k − ∇2η Êk − D(Dk ),∂ta cos φ ∂λa cos φ ∂φ
1 ∂uk Tk′1∂
∂Tk N G= −−(vk Tk′ cos φ) + Ĥk − D(Tk ),∂ta cos φ ∂λa cos φ ∂φ"#(mη̇)k−1/2 T̂k−1/2 − Tk(mη̇)k+1/2 Tk − T̂k+1/2′
Ĥk = Tk Dk −
+
ps
∆σk
ps
∆σk

(2.84)(2.85)(2.86)+κ̂k Tv,k vk · ∇π
X
X
αk
βk
−
Tv,k
vl · ∇π∆Bl −Tv,kvl · ∇π∆Bl
∆σk
∆σk
K

K

l=k

−

αk ′
T
∆σk v,kÊk = Ek +

K
X

K
X

l=k+1

Dl ∆σl −

l=k

βk ′
T
∆σk v,kK
X

Dl ∆σl +l=k+1

Qk + (Qdif f )k, (2.87)
Cp

Wkl (Tv,l − Tl ),(2.88)

k=1

where the vector and matrix of the gravitational wave term (underlined) are
Ck = ∆σk ,(2.89)

Wkl = Cp αl δk≥l + Cp βl δk−1≥l ,(2.90)

Gk = RT̄ ,T̄hkl =[αk ∆σl δk≥l + βk ∆σl δk+1≤l ] .
∆σk(2.91)

Here, δk≤l is 1 if k ≤ l is valid and 0 otherwise.

25(2.92)

We now use the following expressions for time diﬀerences:
1
X t+∆t − X t−∆t ,
2∆t

1
≡
X t+∆t + X t−∆t
2
= X t−∆t + δt X∆t.

δt X ≡
Xt

Then, applying the semi-implicit method to the system of equations, we get ∂πt−C·D ,
δt π =∂t N G∂Dt− ∇2η (ΦS + W T + Gπ t ) − DM (Dt−∆t + 2∆tδt D),
δt D =∂t N G∂Ttδt T =− hD − DH (Tt−∆t + 2∆tδt T).
∂t N G(2.93)(2.94)(2.95)(2.96)(2.97)Thus,t(1 + 2∆tDH )(1 + 2∆tDM )I − (∆t)2 (W h + (1 + 2∆tDM )GCT )∇2η D


∂D
t−∆t
= (1 + 2∆tDH )(1 + ∆tDM )D+ ∆t∂t N G



 
∂T
2
t−∆t
− ∆t∇η (1 + 2∆tDH )ΦS + W (1 + 2∆tDH )T+ ∆t∂t N G  
∂π.
+(1 + 2∆tDH )G π t−∆t + ∆t∂t N G(2.98)

Since the spherical harmonic expansion is used, we can rewrite ∇2η as the following:∇2η = −

n(n + 1),a2

(2.99)tt+∆t =
which enables us to solve the above equations for Dm
n . Then, using (2.95), (2.97) and D
t
t−∆t
t+∆t
2D − D, we can obtain the value of prognostic variables X̂
at t + ∆t.

2.4.4

Time Scheme Properties and Requiments for Time Steps

Let us consider solving the advection equation with the leap-frog method:
∂X
∂X
=c
.
∂t
∂x

(2.100)Assuming X = X0 exp(ikx), the descretized form of the above equation becomes:
X n+1 = X n−1 + 2ik∆tX n .

(2.101)

Assuming X evolves exponentially, we can deﬁne λ such that
λ = X n+1 /X n = X n /X n−1 ,(2.102)2(2.103)λ= 1 + 2ikc∆tλ .

Deﬁning p ≡ kc∆t, the solution becomes:
λ = −ip ±p
1 − p2 .
26(2.104)

The absolute value of those solutions are1 p|p| ≤ 1|λ| =p ± p2 − 1 |p| > 1(2.105)and in the case of |p| > 1, we get |λ| > 1, and the absolute value of the solution increases
exponentially with time. This indicates that the computation is unstable.
In the case of |p| ≤ 1, however, the calculation is neutral since the value of |λ| = 1. However,there are two solutions to λ, one of which, when set to ∆t → 1, leads to λ → 1, while the otherleads to λ → −1, which indicates an osscilating solution. This mode is called “computationalmode” and is one of the problems of the leap frog method. This mode can be damped byapplying a time ﬁlter described later.
Given the horizontal grid spacing ∆x, the maximum value of k becomesπ
max k =.
(2.106)∆x
Then, the condition |p| = kc∆t ≤ 1 requires∆x.
(2.107)πc
In case of a spectral model, using the Earth’s radius a and the maximum wavenumber N ,the requirement becomesa∆t ≤,(2.108)Nc
which is a condition for the numerical stability.
∆t ≤To guarantee the stability of the integration, one needs to take the time step ∆t smaller
than that required by the fastest-propagating mode. If the semi-implicit scheme is not used,
the propagation speed of gravity waves, which can be as fast as 300 m/s, sets the criterion forstability. With the gravity waves taken account of by the semi-implicit method, however, thefastest mode usually becomes the maximum easterly wind Umax . Therefore,a∆t ≤.
(2.109)
N UmaxIn practice, this is multiplied by a factor smaller than 1 for further safety.
2.4.5

Handling of the Initiation of Time Integration

When starting from an initial condition that is not calculated by MIROC 6.0 itself, it is notpossible to give values of all prognostic variables at two time steps t and t − ∆t consistently
with the model dynamics. However, giving inconsistent values for t − ∆t results in a large
computation mode.
To avoid this, a special procedure is followed at the initiation of time integration. Firstly,
assuming X ∆t/4 = X 0 , a 1/4-step integration is performed to obtain X ∆t/2 :X ∆t/2 = X 0 + ∆t/2Ẋ ∆t/4 = X 0 + ∆t/2Ẋ 0 .

(2.110)Then, a 1/2-step integration is performed to yield X ∆t :X ∆t = X 0 + ∆tẊ ∆t/2 .

(2.111)X 2∆t = X 0 + 2∆tẊ ∆t .

(2.112)

Finally, in the normal time step,

From here on, the leap-frog method is executed in the usual manner.
27

2.5

Tracer Advection Scheme

2.5.1

Introduction of Tracer Advection Scheme

MIROC6 adopts a spectral method based on Spherical harmonic expansion to dynamic core.
The spectral method is an excellent method, but it has some drawbacks.
1. Because of Gibbs phenomenon, noisy oscillations are produced when representing a nonsmooth ﬁeld.
2. Associated with Gibbs phenomenon, negative value may occure on grids where they are
not supposed to such as, eg., speciﬁc humidity.
3. Global conservation of conservative quantity is good enough, but local conservation does
not always hold.
4. The property that information is transmitted from upstream to downstream is not always
satisﬁed. In spherical model, information travels instantly to the other side of the world.
Despite of these disadvantages, MIROC has adoptted spectral method as dynamic core. Gibbs
phenomenon usually doesn’t cause any problems. However, when describing the transport of materials with strong discontinuity, the noisy oscillation and unexpected negative values sometimes
appear. For example, water vapor in polar region and the stratosphere often shows discontinuous distribution because there is very small amount of water vapor. Tracers such as aerosols arealso distributed locally and often show large discontinuity. These tracers are easy to aﬀectedby Gibbs phenomenon. Therefore, in MIROC6, water vapor transport and tracer transport arecalculated using a ﬂux form semi-Lagrange (FFSL) scheme (Lin and Rood 1996) instead of using
the spectral method.
Merits of this scheme are described below.
1. Gibbs phenomenon doesn’t occur because it’s based on gridpoint method. Therefore,
non-smooth ﬁelds can be represented with better accuracy.
2. Negative values of tracer quantity can be avoided even in unsmooth ﬁelds.
3. No new extreme values are created.
4. Information is transmitted from upstream to downstream.
5. Conservation is satisﬁed locally and globally.
6. Problems which is induced by narrow grid range in polar region can be avoided.
In the next section, the principle of the tracer advection scheme is introduced in detail, and in
the following section, we describe the actual implementation of the tracer advection scheme.
2.5.2

Principle of the Tracer Advection Scheme

2.5.2.1

Transport Equation in Flux Form

The winds and the tracer distributions are staggered in the Arakawa C-grid (Mesinger andArakawa 1976). As example, three dimension advection equation in (x,y,p) rectangukar coordinate system is given as follows.
∂q
∂q
∂q
∂q
= −u
−v
−ω
∂t
∂x
∂y
∂p
28

(2.113)

Here, q is the amount of tracer (ex. speciﬁc humidity for water vapor), u, v is zonal and meridional velocity respectively. By substituting continuity equation to this, we get the advection
equation in ﬂux form.
∂q
∂
∂
∂
∂
∂ y
∂ p
= − (uq) −(vq) −(ωq) = − F x −
F −
F
∂t
∂x
∂y
∂p
∂x
∂y
∂p

(2.114)Discretizing by x = xi (i = 1, 2, 3...), y = yj (j = 1, 2, 3...), p = pk (k = 1, 2, 3...), the advection
equation is rewritten as∂qi,j,k
1
1
1
y
p
y
p
x
x
=
(Fi−
(Fi,j−(Fi,j,k−− Fi+)+− Fi,j+
1
1
1
1 − F
1 ) +
1)
,j,k
,j,k
,k
,k
i,j,k+
22∂t∆xi,j,k∆yi,j,k∆pi,j,k2
2
2
2(2.115)xHere, Fi−is the ﬂux in x direction at boundary between (i, j, k) and (i − 1, j, k), ∆xi,j,k
1
,j,k
2

is x-direction width of grid represented as (i, j, k). This ﬂux formed equation automaticallyxsatisﬁes the conservation law. The accuracy of the scheme depends on how Fi−
is chosen.
1
,j,k
2

x
Next, how Fi−
is determind in MIRIOC6 is explianed in one dimension in x-direction, for1,j,ksimplicity.
2.5.3

2

The Piecewise Parobolic Method Scheme

In semi-Lagrange scheme, the ﬂux at point xi+ 1 at time t is calculated by using q at point2
xi+ 1 − u∆t at time t − ∆t. The velocity u = u(t) is used assuming constant locally during this2

period, and q at time t − ∆t are used. If CFL condition (| u∆t∆x | < 1) is satisﬁed and ui+ 21 > 0,
xi+ 1 − u∆t is at a point inside grid i.
2
As the value of q at point xi+ 1 − u∆t, qi , which is the average value of point i, can be chosen,
2
assuming that q is constant in the grid.
However, the value of q shows large discontinuity at i + 21 , which is the boundary between
grid i and i + 1 in this assumption. This discontinuity strengthens numerical viscosity, and isunwanted for numerical experiments. Therefore, we want to give some kind of distribution to q,
which is assumed to be constant in a grid, in order to eliminate the discontinuity and enable us
to calculate q at xi+ 1 − u∆t by interpolation. Given distribution must be satisﬁed a condition2
as follows.
Z x 1
i+ 2
1
qi =
q(x)dx(2.116)∆xi x 1i− 2The old editions of MIROC adoptted van Leer method, in which interpolation function is alinear function, but MIROC6 adopts The Piecewise Parobolic Method (PPM) scheme (Colella
and Woodward 1984) , in which interpolation function is a quadrative function (Fig.2). The
FFSL scheme which adopts PPM scheme is called FFSL-3 (Lin and Rood 1996).
In PPM scheme, the distribution is determined as follows.
q(x) = qL,i + ξ(∆qi + q6,i (1 − ξ))x − xi− 12ξ=, xi− 1 ≤ x ≤ xi+ 122∆xi(2.117)Here, qL,i is deﬁned as limx→xi+ 1 = qL,i . qR,i is deﬁned as limx→xi+ 1 = qR,i as well. In PPM
2

2

scheme, q is continuous at boundary i + 12 , therefore qL,i+1 = qR,i = qi+ 1 is hold. In addition,2∆qi = qR,i − qL,i ,1q6,i = 6(qi − (qL,i ) + qR,i )2
29(2.118)Figure 2: The image of interpolation function in The piecewise parobolic method (PPM) scheme.
The interpolation function is in the solid line, the grid mean value is in the dot line.
For calculationg qi+ 1 by interpolation, a ﬁnite integration of q given as follows is introduced.
2
Z x
A(x) =
q(x′ )dx′(2.119)At boundary of grid,A(xi+ 1 ) = Ai+ 1 =2

X2

qk ∆xk(2.120)k≤i

qi+ 1 is calculated by discretization of qi+ 1 = dA/dx|xi+ 1 by using (Aj+k+ 1 , xj+k+ 1 ), k =
2

2

2

2

0, ±1, ±2. Speciﬁcally, qi+ 1 is calculated as follows.

22

qi+ 1 =qi + ∆xi
2

qi+1 − qi1+ Pi+2∆xi+1 + ∆xi
k=i−1 ∆xk

h 2∆x ∆x∆xi + ∆xi−1∆xi+2 + ∆xi+1−)(qi+1 − qi )∆xi+1 + ∆xi ∆xi+1 + 2∆xi
2∆xi+1 + ∆xi∆xi + ∆xi−1∆xi+2 + ∆xi+1 i− ∆xi
δqi+1 + ∆xi+1
δqi∆xi+1 + 2∆xi
2∆xi+1 + ∆xi
×

i

i+1

(

In case the grid width is equal in all grids, Eq.(2.121) can be simply rewritten as1
1
qi+ 1 = (qi−1 + qi) − (δqi − δqi−1 )226Here, δqi is given as∆xi
2∆xi−1 + ∆xi∆xi + 2∆xi+1
δqi =(qi+1 − qi ) +(qi − qi−1 )∆xi−1 + ∆xi + ∆xi+1 ∆xi+1 + ∆xi∆xi−1 + ∆xi(2.121)(2.122)(2.123)

However, in this case, the interpolation function may have extremes in the grid and may notsatisfy monotonicity. In order to avoid such a situation, qi+ 1 should be between qi as qi+1 , and2
δqi is modiﬁed as follows for that.
δm qi = min(|δqi |, 2|qi − qi−1 |, |qi+1 − qi |)if=0(qi+1 − qi )(qi − qi−1 ) > 0,otherwise(2.124)This δqi is used in Eq.(2.121) to calculate qi+ 1 .
2When q(x) is interpolated as Eq.(2.117), by using Courant number deﬁned asC=ui+ 1 ∆t2

(2.125)∆xi+1

x
ﬂux Fi+
1 is wriiten as follows.
2(
ui+ 1 [qR,i − C2 (∆qi − (1 − 32 C)q6,i )](ui+ 1 ≥ 0)x
2
2
Fi+ 1 =
C2
2
ui+ 1 [qL,i+1 + 2 (∆qi+1 + (1 − 3 C)q6,i+1 )] (ui+ 1 ≤ 0)
2

2

30

(2.126)

2.5.3.1

Devices for Taking Long Time Steps

The above argument is stable only if C < 1 When the grid method is adapted to sphericalcoordinate, ∆x is very small in polar region. Therefore, we have to take very small ∆t to satisfy
CFL condition. The special treatment when C > 1 for taking larger ∆t are described below.
Although this method is used for any grid widths, in this subsection we assume that ∆x does
not depend on i for simplicity.
Courant number can be divided into a integral fraction and a decimal fraction.
C = IC + Ĉ,−0.5 ≤ Ĉ ≤ 0.5(2.127)∆xi∆t(2.128)∆xi∆t(2.129)

IC : an integer fraction,

When IC > 0
x
x ˆ
Fi−
+
1 = F
i−IC + 21
2

i
X

qi′

i′ =i+1−IC

When IC < 0i+|IC |
x
x ˆ
Fi−
+
1 = F
i+|I |+ 1
C

2

2

X

i′ =i+1qi′

x ˆ
Where Fi−I
is the ﬂux of point (i − IC + 12 ) calculated by using Ĉ.
+1
C2

As indicated above, in the case the ﬂuid moves in multiple grids during ∆t, we can avoid
instability of numerical calculation by evaluating the ﬂux using the quantity qi′ corresponding
to each grids passed. In actual, these argument is applied only to zonal ﬂux, which can break
CFL condition.
2.5.3.2 The Treatment of Cross Terms
In the case velocity of the ﬂuid is not only in the x-direction or y-direction, only adding
the ﬂux contributions in the x- and y-directions together underestimate the eﬀect of diagonal
advection. To take these cross term into considering, the following procedure is taken. Here, we
discuss this in two-dimensional space, not in one-dimensional.
xWhen calculating x-direction ﬂux Fi+
1 , upstream value of q in y-direction is used as value,j2

of q. That is expressed by the following equation.
1yqi,j= q(xi , yi − vi,j ∆t) + qi,j2

(2.130)Here, q(xi , yi − vi,j ∆t) is calculated by linear interpolation of the two nearest grid points. In thexsame way, when calculating y-direction ﬂux Fi+
1 ,
,j
2

1
x
qi,j
= q(xi − ui,j ∆t, yi ) + qi,j2

(2.131)is used as q.
In the case of three dimensional tracer advection, this procedure is conducted in two dimension.
2.5.4

Actual Tracer Advection Scheme in MIROC6

In this subsection, the procedure implemented in MIROC6 of the tracer advection scheme isdescribed. Although MIROC6 adopts σ-p hybrid coordinate as vertical coordinate, the tracer
advection scheme is largely based on σ coordinate because previous version of MIROC adopted
σ coordinate. Therefore, ﬁrstly the procedure under σ coordinate system is described. After
this, the changes in the hybrid coordinate system from the σ coordinate system is described.
31

2.5.4.1 σ-Coordinate
The transport equation in σ coordinate on the sphere is expressed as∂P S q
∂t

1
∂
1
∂
∂
(P S uq) −(P S vq cos φ) −(P S σ̇q)a cos φ ∂λa cos φ ∂φ
∂σ
1
∂
1
∂
∂
(F λ ) −(F φ ) −(F σ )
a cos φ ∂λa cos φ ∂φ
∂σ

= −
=

(2.132)P S is surface pressure, q is quantity of tracers. Continuity equation is given by considering thecase of q = 1.
∂P S
1
∂
1
∂
∂
=−
(P S u) −(P S v cos φ) −(P S σ̇)∂ta cos φ ∂λa cos φ ∂φ∂σ(2.133)

Assuming that grid is equally spaced in zonal direction, the transport equation is discretized asfollows.
S q∂P,i,j,k
i,j,k

∂t

=

1
−Gφ
))+(Gσi,j,k− 1 −Gσi,j,k+ 1 )] (2.134)[(Gλi− 1 ,j,k −Gλi+ 1 ,j,k )+(Gφi,j+ 21 ,ki,j− 12 ,k
2
2
2
2
∆Dj,k

Here,
λ
Gλi− 1 ,j,k = Fi−∆yj ∆σk = (P S uq)i− 1 ,j,k ∆yj ∆σk1,j,k(2.135)φs
Gφ= Fi,j−1 ∆xj− 1 ∆σk = (P vq)i,j− 1 ,k ∆xj− 1 ∆σki,j− 1 ,k,k(2.136)ηS
Gσi,j,k− 1 = Fi,j,k−
1 ∆xj ∆yj = (P σ̇q)i,j,k− 1 ∆xj ∆yj(2.137)

2

2

2

2

2

2

2

2

2

2

2

And
∆Dj,k = a cos φj ∆λ∆φj ∆σk ,∆xj = a cos φj ∆λ,∆yj = a∆φj

(2.138)This ﬂux form equation ensure the conservation.
For the calculation of the time-averaged mass ﬂux across the cell boundary, the winds andthe tracer distributions are staggered in the Arakawa C-grid (Mesinger and Arakawa 1976). The
horizontal winds at the cell boundary, ui− 1 ,j,k , vi− 1 ,j,k , are reconstructed by using the mass2
2
convergence ﬁeld in the spectral model and the discretized continuity eqution:
S∂Pi,j,k

=

∂t

1
φ
φ
σ
λ
σ
[(V λ 1− Vi+
) + (Vi,j−− Vi,j+11 )) + (Vi,j,k− 1 − Vi,j,k+ 1 )]
1
,j,k
,k
,k
2
2
2
∆Dj,k i− 2 ,j,k2
2(2.139)λ,V φdenote zonal, meridional, and vertical mass-weighted wind at theHere, Vi−,V σ
1,j,k i,j− 1 ,k i,j,k− 1222cell boundary, rspectively. That is,λVi−= (P S u)i− 1 ,j,k ∆yj ∆σk1,j,k(2.140)φ
Vi,j−= (P S v)i,j− 1 ,k ∆xj− 1 ∆ηk
1,k(2.141)σSVi,j,k−1 = (P σ̇)i,j,k− 1 ∆xj ∆yj(2.142)

2

2

2

2

2

2

2

∆Dj,k denotes the cell volume, and ∆xj , ∆yj , and ∆σk denote zonal, meridional and vertical
width of the cell, respectively. That is ∆Dj,k = a cos φj ∆λ∆φj ∆σ, ∆xj = a cos φj ∆λ and∆yj = a∆φj .
The following are the procedure for the calculation of tracer advection in the staggeringgrided horizontal and vertical wind ﬁelds:32

1. Surface pressure P S (t + ∆t) and horizontal wind v(t + ∆t) are predicted in the spectral
model.
2. The horizontal component of mass ﬂux divergence at time step t is calculated by using
spherical harmonics. The mass ﬂuxes at time step t are reconstructed from the values at
t + ∆t and t − ∆t because MIROC applies semi-implicit scheme for the time-integrationof surface pressure. Zonal and meridional component of mass ﬂux divergence are:Cx = −

1∂(P s u),a cos φ ∂λCy = −

1∂(P s v cos φ)a cos φ ∂λ(2.143)

3. By using Cx and Cy , V λ , V φ , V σ are calculated as follows.
λ
λ
x
Vi−
− Vi+
= Ci,j,k∆Dj,k ,
1
1
,j,k
,j,k
2

2

y
λ
λ
Vi,j−
− Vi,j+= Ci,j,k
∆Dj,k
1
1
,k
,k
2

(2.144)

2

The boundry conditions are V φ = 0 at the North Pole and South Pole, σ = 1 at surface
and V σ = 0 at σ = 0. The condition for V λ = 0 is:
X
X
λ
S
Vi−
=
Pi,j,kui,j,k ∆yj ∆σk(2.145)
1
,j,k
i

2

i

That means zonal mean of zonal mass transport is equal to that in the spectral modelgrid. Here, the following equation must be satisﬁed for boundary condition V φ = 0 at theNorth Pole and the South Pole.
X y(2.146)
Ci,j,k ∆Dj,k = 0
jHowever, this is not always satisﬁed (On the other hand,within numerical error.).

P P
i

y
j Ci,j,k ∆Dj,k = 0 is valid

In order to satisfy the boundary condition, the following correction is made.
y
y
− δC,
← Ci,j,kCi,j,kHere, δC =

y
x
+ δC
← Ci,j,kCi,j,k

(2.147)PPyηj Ci,j,k ∆Dj,k /j ∆Dj,k . Vertical velocity V is obtained by using
S
X∂Pi,j,k∂t∆Dj,k =kX

yx
+ Ci,j,k)(Ci,j,k(2.148)k(The contents so far are in [TRACEG] of dtrcr.F. The rest of the content is in [GTRACE]of dtrcr.F.)
4. Gλ , Gφ , Gσ are calculated by PPM scheme from V λ , V φ , V σ .
s q
5. Pi,j,ki,j,k at time step t + ∆t is calulated by integration of Eq.(2.134) by leap frog methodfrom Gλ , Gφ , Gσ .
s
6. qt+∆t is calculated by dividing (P s q)t+∆t by Pt+∆t. There is small quantity of diﬀerencessbetween Pt+∆t from Eq. (2.148) and Pt+∆t in the spectral model, because semi-implicitstime integration scheme is applied. Pt+∆tfrom Eq. (2.148) is applied at present forthe consistency of mass advection. Mass Conservation is not strictly satisﬁed because ofs
the discrepancy between the surface pressure in the spectral model and from Pt+∆t
Eq.
(2.148).

33

2.5.4.2 σ-p Hybrid Coordinate
The transport equation in η coordinate (σ − p hybrid coordinate) on the sphere is:
∂mq
∂t

1
∂
1
∂
∂
(muq) −
(mvq cos φ) −(mη̇q)
a cos φ ∂λa cos φ ∂φ
∂η
1
∂
1
∂
∂
(F λ ) −(F φ ) −(F η )
a cos φ ∂λa cos φ ∂φ
∂η

= −
=

(2.149)∂pHere, m corresponds to the density of the coordinate and is deﬁned as m = ∂η
. if you lookat Eq. (2.132), you can ﬁnd that diﬀerence of σ coordinate and η coordinate is only that P S
replaces m. The actual tracer advection in η coordinate is mostly the same as σ coordinate.
In the scheme in η coordinate, the following variables are calculated in the same way as theway Gλ , Gφ , Gη is calculated in σ coordinate, except ∆σ replaces with ∆η and σ̇ replaces withη̇.

G′λ= (P S uq)i− 1 ,j,k ∆yj ∆ηk ,i− 1 ,j,k2

2G′φ= (P S vq)i,j− 1 ,k ∆xj− 1 ∆ηk ,i,j− 1 ,k2

2G′η= (P S η̇q)i,j,k− 1 ∆xj ∆yj
i,j,k− 12

2

2(2.150)S
λφηIn the time integration step, multiplying G′ by m/P , G , G , G is calculated. After that, mq
at time step t + ∆t is calculated by leap-frog method as well as σ coordinate.
In actual source code, conbining to dividing by m to calculate q at time step t + ∆t, q at
point (i, j, k) in time step t + ∆t is calculated as follows.
q t+∆t =S,t−∆t∆Ak + ∆Bk Pi,j,k2∆tq t−∆t +S,t+∆t i,j,k∆D∆Ak + ∆Bk Pi,j,k× [(G′λ,t− G′λ,t) + (G′φ,t− G′φ,t)) + (G′η,t− G′η,t)]i− 1 ,j,k
i+ 1 ,j,ki,j− 1 ,ki,j+ 1 ,ki,j,k− 1
i,j,k+ 1
2

×

2

2

2

S,t
∆Ak + ∆Bk Pi,j,k1

S,tPi,j,kS,t+∆t∆Ak + ∆Bk Pi,j,k2(2.151)2

Here,A, B is the coeﬃcients for η coordinate, ηk+ 1 = Ak+ 1 /p0 + Bk+ 1 and ∆Ak = Ak− 1 −
2

2

2

2

S
= ∆pi,j,k (More details in the section of
Ak+ 1 , ∆Bk = Bk− 1 − Bk+ 1 . And ∆Ak + ∆Bk Pi,j,k2
22
the vertical discretization).
2.5.4.3 The Mass Fluxes into/out of Polar Caps

The mass ﬂuxes into/out of polar caps are calculated by using the semi–Lagrangian scheme
in the polar stereo projection (cf. Fig.3). The horizontal average at the highest latitude band isassumed to be preserved before/after ﬂux calculation for the mass conservation. The sequenceof calculation is:
1. Zonal average of P S q at time step t is calculated at the highest latitude band (j = jN , jS ),and is assumed to equal P S q at the pole.
2. Horizontal wind at the highest latitude bands is projected into the orthogonal coordinate
system centering around the pole, and q at time step t + ∆t is estimated by using the value
at the departure point .
3. Zonal anerage of P S q at time step t + ∆t is ﬁxed to that at t.

34

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

35

2.6

Summary of the Dynamical Core

In this section, we enumerate the calculations performed in the dynamical core, although they
overlap with the previous descriptions.
2.6.1

Conversion of Horizontal Wind to Vorticity and DivergenceObtain grid point values of vorticity and divergence from the grid point values of uij , vij forhorizontal wind. First, we obtain the vorticity and divergence in spectral space, ζnm , Dnm ,
1
ζnm =
I

J
I X
X
i=1 j=1

wj
wj
1 XX
1 XX
∂
imuij cos φj Ynm∗ ij−
vij cos φj (1 − µ2j ) Ynm∗ ij(2.153);.
2I∂µa(1 − µj ) Ia(1 − µ2j )
I

Dnm =

J
I
wj
wj
1 XX
∂
+
(2.152),uij cos φj (1 − µ2j ) Ynm∗ ij
2∂µa(1 − µj ) I i=1 j=1a(1 − µ2j )imvij cos φj Ynm∗ ij

JI

i=1 j=1J

i=1 j=1

The grid point value is calculated by

ζij = ReNX

NX

ζnm Ynm ij ,(2.154)m=−N n=|m|and so on.
Corresponding ﬁle & subroutines: [G2Wpush, G2Wtrans, G2Wshift, W2Gpush, W2Gtrans,
W2Gshift (xdsphe.F)]
2.6.2

Calculating a Virtual Temperature

virtual Temperature Tv is ,Tv = T (1 + ϵv q − l) ,(2.155)

However, it is ϵv = Rv /R − 1 and Rv is the gas constant for water vapor (461 Jkg−1 K−1 )and R is the gas constant for air (287.04 Jkg−1 K−1 ).
Corresponding ﬁle & subroutine: [VIRTMD (dvtmp.F)]
2.6.3

Calculating the Pressure Gradient Term

The pressure gradient term ∇π = p1S ∇pS is ﬁrst used to deﬁne the πnm
1 XX(ln pS )ij Ynm∗ ij wj ,
I
I

πnm =

J

(2.156)i=1 j=1

to a spectral representation and then ,1a cos φ



∂π
∂λ



N
X
1
=
Re
a cos φij

NX

m=−N n=|m|36

imX̃nm Ynm ij ,(2.157)

1
a



∂π
∂φ



N
X
1
Re
a cos φ=ij

NX

πnm (1 − µ2 )m=−N n=|m|∂ m
Y.
∂µ n ij(2.158)Corresponding ﬁle & subroutine: [PSDOT (dgdyn.F)]
2.6.4

Diagnosis of Vertical Dlow

Pressure change term, and lead DC,X∂π=−{Dk ∆σk + (vk · ∇π)∆Bk }∂t(2.159)K(mη̇)k−1/2∂π X= −Bk−1/2−{Dl ∆σl + (vl · ∇π)∆Bl }ps∂t(2.160)

K

k=1

l=k

and its non-gravity components.
G(mη̇)N
k−1/2

ps

∂π
∂t

N G
=−

K
X

vk · ∇π∆Bk(2.161)k=1

= −Bk−1/2

∂π
∂t

N G
−

K
X

vl · ∇π∆Bl(2.162)l=kCorresponding ﬁle and subroutine: [PSDOT (dgdyn.F)]
2.6.5

Tendency Terms due to Advection

Momentum advection term:(mη̇)k+1/2 uk − uk+1(mη̇)k−1/2 uk−1 − uk(Au )k =(ζk + f )vk −+ps∆σk−1 + ∆σkps∆σk + ∆σk+1
1 ∂π−(Cp Tv,k κ̂ − RT̄ ) + Fx
(2.163)a cos φ ∂λ(mη̇)k−1/2 vk−1 − vk(mη̇)k+1/2 vk − vk+1
(Av )k =−(ζk + f )uk −+ps∆σk−1 + ∆σkps∆σk + ∆σk+1
1 ∂π−(Cp Tv,k κ̂ − RT̄ ) + Fy
(2.164)a ∂φTemperature advection term:(uT ′ )k = uk (Tk − T̄ )37

(2.165)

(vT ′ )k = vk (Tk − T̄ )(2.166)"Hk =(mη̇)k−1/2 T̂k−1/2 − Tk(mη̇)k+1/2 Tk − T̂k+1/2
Tk′ Dk −
+
ps
∆σl
ps
∆σl
+ κ̂k vk · ∇πTv,k − αkKX(Dl ∆σl + (vl · ∇π)∆Bl )

l=k

− βk

K
X

(Dl ∆σl + (vl · ∇π)∆Bl )

l=k+1

#

Tv,k
∆σk

Tv,k
∆σk

(2.167)

Water vapor advection term:

Rk = qk Dk −(uq)k = uk qk(2.168)(vq)k = vk qk(2.169)(mη̇)k+1/2 qk − qk+1
1 (mη̇)k−1/2 qk−1 − qk
+
2
ps
∆σk
ps
∆σk

(2.170)

Corresponding ﬁle & subroutine: [GRTADV, GRUADV (dgdyn.F)]
2.6.6

Transformation of Prognostic Variables to Spectral Space(122) and (123).
t−∆tTransform ut−∆t, vij
to a spectral representation of vorticity and divergence ζnm , Dnm . Furij
thermore, transforming the temperature T t−∆t , speciﬁc humidity q t−∆t , and π = ln pt−∆t
to
S

1 XX
Xij Ynm∗ ij wj ,
I
I

Xnm =

J

(2.171)i=1 j=1

to a spectral representation.
Corresponding ﬁle & subroutine: [G2Wpush, G2Wtrans, G2Wshift (xdsphe.F)]
2.6.7

Transformation of Tendency Terms to Spectral Space

Tendency Term of Vorticity
wj
1 XX
im(Av )ij cos φj Ynm∗ ijIa(1 − µ2j )
I

∂ζnm
=
∂t

J

i=1 j=1

wj
1 XX∂(Au )ij cos φj (1 − µ2j ) Ynm∗ ijI∂µa(1 − µ2j )I+

J

i=1 j=138(2.172)

The non-gravity wave component of the tendency term of the divergence


∂Dnm
∂t

N G

wj
1 XX
im(Au )ij cos φj Ynm∗ ijIa(1 − µ2j )I=

J

i=1 j=1

wj
1 XX∂(Av )ij cos φj (1 − µ2j ) Ynm∗ ijI∂µa(1 − µ2j )I−

J

i=1 j=1n(n + 1) 1 X X
Êij Ynm∗ ij wj
a2
I
I

−

J

(2.173)i=1 j=1

The non-gravity wave component of the tendency term of temperature


∂Tnm
∂t

N G

wj
1 XX
im(uT ′ )ij cos φj Ynm∗ ijIa(1 − µ2j )
I

−

=

J

i=1 j=1

wj∂
1 XX(vT ′ )ij cos φj (1 − µ2j ) Ynm∗ ijI∂µa(1 − µ2j )I+

J

i=1 j=1

1 XX
Ĥij Ynm∗ ij wj
I
I

+

J

(2.174)i=1 j=1

Tendency term of water vaporwj
1 XX
im(uq)ij cos φj Ynm∗ ijIa(1 − µ2j )
I

∂qnm
=
∂t

−

J

i=1 j=1

+

1
I

I X
J
X

(vq)ij cos φj (1 − µ2j )i=1 j=1

1 XX
Rij Ynm∗ ij wj
I
I

+

wj
∂ m∗
Yn ij
∂µ
a(1 − µ2j )J(2.175)i=1 j=1

Corresponding ﬁle & subroutines: [G2Wpush, G2Wtrans, G2Wshift (xdsphe.F)]
2.6.8

Time Integration in Spectral Space

Equations in matrix formt

(1 + 2∆tDH )(1 + 2∆tDM )I − (∆t)2 (W h + (1 + 2∆tDM )GCT )∇2σ D


∂D
t−∆t
= (1 + 2∆tDH )(1 − ∆tDM )D+ ∆t∂t N G ∂T−∆t∇2σ (1 + 2∆tDH )ΦS + W (1 − 2∆tDH )Tt−∆t + ∆t∂tG  N
∂π+(1 + 2∆tDH )G π t−∆t + ∆t.
∂t N G39(2.176)

Using LU decomposition, D̄ is obtained by solving for
∂T
=
∂t
∂π
=
∂t





∂T
∂t

∂π
∂t


− hD

(2.177)−C·D(2.178)NGNG

Calculate the value of the spectrum in ∂T/∂t, ∂π/∂t and then calculate the value of thespectrum in t + ∆t using

ζ

t+∆t



∂ζ
t−∆t
= ζ
+ 2∆t
(1 + 2∆tDM )−1∂t

Dt+∆t = 2D̄ − Dt−∆t∂TT t+∆t = T t−∆t + 2∆t(1 + 2∆tDH )−1
∂t


∂q
t+∆t
t−∆t
(1 + 2∆tDE )−1
q
= q
+ 2∆t
∂t
∂π
π t+∆t = π t−∆t + 2∆t∂t(2.179)(2.180)(2.181)(2.182)

Corresponding ﬁle & subroutine: [TINTGR (dintg.F)]
2.6.9

Transformation of Prognostic Variables to Grid Point Values

Obtain grid values of horizontal wind speed from the spectral values of vorticity and divergence(ζnm , Dnm ) uij , vij .

N
N 
X
X
1
ima
a
m
2 ∂
m
m m
uij =
ζ (1 − µ ) Yn ij −D YRecos φjn(n + 1) n∂µn(n + 1) n n ij

(2.183)


N
N 
X
X
1
ima
a
m m
m
2 ∂
m
vij =
Re
−
ζ Y
−
D̃ (1 − µ ) Yn ij
cos φjn(n + 1) n n ij n(n + 1) n∂µ(2.184)m=−N n=|m|n̸=0m=−N n=|m|n̸=0Furthermore,

Tij = ReNX

N
X

Tnm Ynm ij ,(2.185)m=−N n=|m|Tij , πij , qij , and so on,pS ij = exp πij
to calculate.
Corresponding ﬁle & subroutines: [W2Gpush, W2Gtrans, W2Gshift (xdsphe.F)]40(2.186)

2.6.10

Diﬀusion Correction along Pressure Level

The horizontal diﬀusion is applied on the surface of η−plane, but it can cause problems in large
slopes, such as transporting water vapor uphill and causing false precipitation at the top of amountain. To mitigate this problem, corrections have been made for T, q, l to make the diﬀusion
closer to that of the p surface, e.g., for T, q, l.
∂T∂σD(−1)ND /2 K∇N
η p·
∂p
∂σ
∂T
ND /2
D
D
= (−1)ND /2 K∇NK∇N
η T − (−1)
η π·σ
∂σ
∂T
= D(T ) − D(π)σ(2.187)
∂σ

ND /2
D
D
Dp (T ) = (−1)ND /2 K∇NK∇N
p T ≃ (−1)
η T −

So,

Tk ← Tk − 2∆tσk

Tk+1 − Tk−1
D(π)σk+1 − σk−1(2.188)and so on. In D(π), the spectral value of π is converted to a grid by multiplying the spectral
value of πnm by the spectral representation of the diﬀusion coeﬃcient.
Corresponding ﬁle & subroutine: [CORDIF (ddifc.F)]
2.6.11

Frictional Heat Associated with Diﬀusion

Frictional heat from diﬀusion is ,QDIF = − (uij D(u)ij + vij D(v)ij )(2.189)It is estimated that Therefore,

Tk ← Tk −2∆t(uij D(u)ij + vij D(v)ij )Cp(2.190)

Corresponding ﬁle & subroutine: [CORDIF (ddifc.F)]
2.6.12

Horizontal Diﬀusion and Rayleigh Friction

The coeﬃcients of horizontal diﬀusion can be expressed spectrally,"DM mn = KM#
2 ND /2
−
+ KR
a2


n(n + 1) ND /2m
DH n = KMa2n(n + 1) ND /2m
DE n = KEa2n(n + 1)a2

ND /2

KR is the Rayleigh coeﬃcient of friction. The Rayleigh coeﬃcient of friction is41(2.191)0KR = KRz − zR
1 + tanhHR(2.192)

However, the proﬁle is given in the same way as However,z = −H ln σ(2.193)0 = (30day)−1 and z = −H ln σ . The standard
The results are approximate to those of KR
topR−1
0 = (30day) , z = −H ln σvalues are KR
top (σtop : top level of the model), H = 8000 m, andR
HR = 7000 m.
Corresponding ﬁle & subroutine [DSETDF (dsetd.F)]

2.6.13Time Filter

To reduce numerical mode associated with leap frog scheme, time ﬁlter is applied every time
step. MIORC6 used modiﬁed Asselin time ﬁlter (Williams, 2009), which is updated versionof Asselin(1972) used previous version of MIROC. Although Asselin time ﬁlter attenuate highfrequency physical mode, bringing low accuracy of leap frog scheme, current time ﬁlter succeededin suppressing it.
Modiﬁed Asselin ﬁlter is expressed as following equation¯ t = X̄ t + να[X̄¯ t−∆t − 2X̄ t + X t+∆t ]X̄

(2.194)¯ t−∆t − 2X̄ t + X t+∆t ]X̄ t+∆t = X t+∆t + ν(1 − α)[X̄(2.195)where bar indicates time ﬁlter. The parameters set to ν = 0.05, α = 0.5. Assuming α = 1,
modiﬁed Asselin ﬁlter is same as Asselin ﬁlter.
In the model,¯ t∗ = (1 − να)−1 [(1 − 2να)X̄ t + ναX̄¯ t−∆t ]X̄

(2.196)

is ﬁrstly calculated at MODULE: [DADVNC] where transformation of prognostic variableto grid
point values. And then, X t−∆t − 2X t is stored. When the X t+∆t is obtained later, time ﬁlter
conduct at MODULE [TFILT],¯ t = (1 − να)X̄¯ t∗ + ναX t+∆tX̄¯ t−∆t − 2X̄ t + X t+∆t ]X̄ t+∆t = X t+∆t + ν(1 − α)[X̄(2.197)(2.198)Corresponding ﬁle & subroutine: [DADVNC (dadvn.F)]
2.6.14

Correction for Conservation of Mass

In the spectral method, the global integral of π = ln pS is preserved with rounding errors
removed, but the preservation of the mass, i.e. the global integral of pS is not guaranteed.
Moreover, a wavenumber break in the spectra sometimes results in negative values of the water
vapor grid points. For this reason, we perform a correction to preserve the masses of dry air,water vapor, and cloud water, and to remove the regions with negative water vapor content.
42

Before entering dynamical calculations, [FIXMAS], the global integrals of water vapor andcloud water are calculated for Mq , Ml .
Mq0 =

X

qpS ∆λi wj ∆σk(2.199)lpS ∆λi wj ∆σk(2.200)

ijk

Ml0 =

X
ijk

In the ﬁrst step of the calculation, the dry mass Md is calculated and stored.
Md0 =X(1 − q − l)pS ∆λi wj ∆σk(2.201)

ijk

After exiting dynamical calculation, [MASFIX], the following procedure is followed.
First, negative water vapor is removed by dividing the water vapor from the grid points
immediately below the grid points. Suppose that qk < 0 is used,qk′ = 0∆pk′
= qk−1 +
qk−1
qk
∆pk−1

(2.202)(2.203)′≥ 0.
However, this should only be done if it is qk−1
Next, set the value to zero for the grid points not removed by the above procedure.

3. calculate the global integral value of Mq and multiply the global water vapor content by
a ﬁxed percentage so that it is the same as that of Mq0 .

q ′′ =

Mq0 ′q
Mq(2.204)

4. correct for dry air mass Likewise calculate Md ,p′′S =

Md0
pS
Md

Corresponding ﬁle & subroutine: [FIXMAS, MASFIX (dmfix.F)]

43(2.205)

2.7

Computational Flow of Dynamical Core

In this section, calculations of dynamical component based on coding are summarized. [module
name(file name)]
2.7.1

Overview of Dynamical Core

1. Calculate dynamical term [DYNTRM(dterm.F)]
1.1 Calculate volticity and divergence on wave space and get grid value. [G2W, W2G(xdsphe.F)]
1.2 Diagnose stream function and potential velocity [DYNTRM(dterm.F)]
1.3 Diagnose surface pressure advection, its tendency & vertical ﬂow [PSDOT(dgdyn.F)]
1.4 Calculate factor for hydrostatic eq. & interporation of temprature on Hybrid coord.
[CFACT(dcfct.F)]
1.5 Calculate virtual temprature [VIRTMD(dvtmp.F)]
1.6 Calculate temperature advection [GRTADV(dgdyn.F)]
1.7 Calculate momentum advection [GRUADV(dgdyn.F)]
1.8 Spectral transform of tendency terms [G2W(xdsphe.F)]
2. Time integration of equation DYNSTP(dstep.F)
2.1 Calculate tracer transport [TRACEG(dtrcr.F)]
2.2 Time integration on wave space [TINTGR(dintg.F)]
2.3 Time integration of tracer terms [GTRACE(dtrcr.F)]
2.4 Time ﬁlter [DADVNC(dadvn.F)]
2.5 Get horizontal wind of grid value from wave space [W2G(xdsphe.F)]
2.6 Correction of pressure-level diﬀusion [CORDIF(ddifc.F)]
2.7 Correction of horizontal friction heating [CORFRC(ddifc.F)]

44

3

Physics

3.1

Overview of Physical Parameterizations

As a physical process, we can consider the following• Cumulus convection• Shallow convection• Large scale condensation• Cloud microphysics• Radiation• Turbulence• Surface ﬂuxes• Gravity wave dragWe compute the time-varying terms Fx , Fy , Q, M, S for the prognostic variables from theseprocesses, and perform time integration.
3.1.1

Time Integration of Physical Parameterizations

NOTE: the descriptions in this section are outdated.
In terms of time integration of predictors, we can classify the physical Parameterizations in
the following three orders of execution.
1. Cumulus convection, shallow convection, large-scale condensation, and cloud microphysics
2. Radiation, turbulence and surface ﬂuxes
3. Gravitational wave dragFor cumulus convection, shallow convection, large-scale condensation, and cloud microphysics, the values are updated by the usual Euler diﬀerence as follows.
T̂ t+∆t,(1) = T̂ t+∆t + 2∆tQCU M (T̂ t+∆t )(3.1)

T̂ t+∆t,(2) = T̂ t+∆t,(1) + 2∆tQLSC (T̂ t+∆t,(1) )(3.2)

Note that the large-scale condensation scheme is updated by the cumulus convection scheme.
In practice, the routines of cumulus convection and large-scale condensation output the heating
rates and so on, and the time integration is performed immediately afterwards by MODULE:[GDINTG].
The calculations of the radiative, vertical diﬀusion, ground boundary layer and surface processes in the following groups are basically performed with these updated values (T̂ t+∆t,(1) , q̂ t+∆t,(2) ,etc.). However, in order to calculate some of the terms as implicit, we calculate the heating rates
and so on for all of these terms together, and then perform time integration at the end. In otherwords, if we write symbolicallyT̂ t+∆t,(3) = T̂ t+∆t,(2) + 2∆tQRAD,DIF,SF C (T̂ t+∆t,(2) , T̂ t+∆t,(3) )(3.3)

That would be the gravitational wave resistance, mass modulation, and dry convection modulation are the same as those for cumulus convection and large-scale condensation.
T̂ t+∆t,(4) = T̂ t+∆t,(3) + 2∆tQADJ (T̂ t+∆t,(3) )
45

(3.4)

3.1.2

Various Physical Quantities

Here are deﬁnitions of various physical quantities that can be computed simply from the prognostic variables. Some of them are calculated with MODULE:[PSETUP].
1. Virtual temperature
Virtual Temperature Tv is calculated as follows.
Tv = T (1 + ϵv q − l)

(3.5)

2. Air density
The air density ρ is calculated as follows.

ρ=

p
RTv

(3.6)

3. Altitude
The altitude z is evaluated in the same way as the calculation of the geopotential height in
the dynamics.

z=

Φ
g

Φ1 = Φs + Cp (σ1−κ − 1)Tv,1
Φk − Φk−1 = Cp

σk−1/2
σk

κ




 
σk−1/2 κ− 1 Tv,k + Cp 1 −Tv,k−1
σk−1

(3.7)

(3.8)

(3.9)

4. Half-level temperature
Half-level temperature is calculated by performing a linear interpolation on ln p, i.e., ln σ.

Tk−1/2 =ln σk−1 − ln σk−1/2
ln σk−1/2 − ln σk
Tk +
Tk−1
ln σk−1 − ln σk
ln σk−1 − ln σk(3.10)

5. Saturated speciﬁc humidity
The saturated speciﬁc humidity q ∗ (T, p) is approximated using the saturated vapor pressuree∗ (T ),q ∗ (T, p) =ϵe∗ (T ).
p

Here, it is ϵ = 0.622,
46(3.11)

L
1 ∂e∗v=∗ev ∂T
Rv T 2(3.12)

Therefore, if the latent heat of evaporation (L) and the gas constant of water vapor (Rv ) are
held constant,Le (T ) = e (T = 273K) exp
Rv
∗

∗



1
1
−
273 T


,

(3.13)

where e∗ (T = 273[K]) = 611[hPa].
From eq.3.12,
L ∗
∂q ∗
=
q (T, p).
∂T
Rv T 2(3.14)

It is noted that if the temperature is lower than the freezing point 273.15[K], the sublimation
latent heat L + LM is used as the latent heat L.
6. Dry static energy and moisture static energy
The dry static energy s is deﬁned by

s = Cp T + gz

(3.15)The moisture static Energy h is deﬁned by

h = Cp T + gz + Lq

47

(3.16)

3.2
3.2.1

Cumulus Scheme
Outline of Cumulus Scheme

The Chikira scheme (Chikira and Sugiyama 2010) has been adopted since version 5 of MIROC.
It represents updrafts, downdrafts, their detrainment and compensating downward motion overthe surrounding area as well as microphysical processes associated with updrafts and downdrafts.
The updraft is based on an entraining plume model, where the mass ﬂux increases upward
due to lateral entrainment. The detrainment occurs only at the cloud top which is deﬁned as theneutral buoyancy level of the updraft air parcel. The lateral entrainment rate is formulated interms of buoyancy and vertical velocity of the air parcel at each level following Gregory (2001).
The momentum transport is formulated following Gregory et al. (1997).
The cloud base mass ﬂuxes are determined by the prognostic convective kinetic energy closure
proposed by Arakawa and Xu (1990) and Xu (1991, 1993), which was adopted in the prognostic
Arakawa–Schubert scheme (Randall and Pan 1993; Pan 1995; Randall et al. 1997; Pan andRandall 1998). The convective kinetic energy increases by buoancy and decreases by dissipation.
The cloud types are spectrally represented according to the updraft vertical velocity at thecloud base. Larger (smaller) vertical velocities give smaller (larger) entrainment rates which
result in higher (lower) cloud tops. The cloud base is diagnosed as the lifting condensation levelof the air parcel at the lowest model layer.
The scheme has a simple downdraft model, where a part of the precipitation caused by the
updrafts evaporates and forms the cold air which enters into the downdrafts. The detrainmentof the downdraft mass ﬂuxes occurs at the neutral buoyancy level and near the surface.
The interaction of the updrafs and downdrafts with the surrounding environment is formulated following Arakawa and Schubert (1974). The areal fractions of the updrafts and downdrafts
are assumed to be suﬃciently small and the grid-mean prognostic variables are supposed to be
the same as those over the environmental area, which are changed by the detrainment of the
updrafts and downdrafts, the compensating subsidence and the evaporation and sublimation of
the precipitation associated with the updrafts.
The input variables to this scheme are temperature T , speciﬁc humidity q, cloud liquid water
ql , cloud ice qi , zonal wind u, meridional wind v, all tracers including aerosols and greenhousegases, height z, pressure p, and cloud cover C. The scheme gives the tendencies of T , qv , ql , qi ,u, v, C and all the tracers. The vertical proﬁles of the rainfall and snowfall ﬂuxes, cloud liquidwater, cloud ice and cloud fraction associated with the updrafts are also output as diagnostic
variables.
The procedure of the calculations is given as follows along with the names of the subroutines.
1. Calculation of cloud base CUMBAS.
2. Calculation of in-cloud properties CUMUP.
3. Calculation of cloud base mass ﬂux CUMBMX.
4. Calculation of cloud mass ﬂux, detrainment, and precipitation CUMFLX.
5. Diagnosis of cloud water and cloud cover by cumulus CUMCLD.
6. Calculation of tendencies by detrainment CLDDET.
7. Calculation of freezing, melting, evaporation, sublimation, and downdraft mass ﬂux CUMDWN.
8. Calculation of tendencies by compensating subsidence CLDSBH.
9. Calculation of cumulus momentum transport CUMCMT.
10. Calculation of tracer updraft CUMUPR.
11. Calculation of tracer downdraft CUMDNR.
12. Calculation of tracer subsidence CUMSBR.
13. Fixing tracer mass CUMFXR .
48

3.2.2

Interaction between Cumulus Ensemble and Large-Scale Environment

Following Arakawa and Schubert (1974), the equations for tendencies of the grid-mean variables
are written as∂ h̄∂ h̄ X=M+
Dj hj (zT,j ) − h̄ ,∂t∂z

(3.17)∂ q̄∂ q̄ X=M+
Dj [qj (zT,j ) − q̄] ,∂t∂z(3.18)jj

where M , D, h denote total mass ﬂux, detrainment mass ﬂux and moist static energy. q
is a substitute for qv , ql and qi and any tracers which are calculated in the same way. zT is
the height of the updraft. The hats indicate in-cloud properties, the overbars grid-mean. The
subscripts j are an index for the updraft types.
The total mass ﬂux M and detrainment D are deﬁned as

M (z) =

X

Mu,j + Md ,(3.19)j

Dj (z) = Mu,j (zT,j )δ(z − zT,j )(3.20)

respectively, where Mu and Md denote mass ﬂuxes of updraft and downdraft respectively.
The updraft mass ﬂux is formulated as

Mu,j (z) = MB,j ηj (z)(3.21)

where MB and η are the updraft mass ﬂux at its cloud base and normalized mass ﬂux.
3.2.3

Cloud Base

The cloud base is determined as the lifting condensation level of the air at the lowest model
layer. It is deﬁned as the smallest z which satisﬁesq¯t (z1 ) ≥ q¯v ∗ +γh̄(z1 ) − h̄∗ (z) ,Lv (1 + γ)

(3.22)

where qt denotes total water, Lv the latent heat of vaporization, z1 the height of the lowest
model layer at the full level and
Lv
γ≡
Cp



∂ q̄ ∗∂ T̄.

(3.23)p̄

Cp denotes the speciﬁc heat of dry air at constant pressure and the stars indicate saturation
values.
The normalized mass ﬂux below the cloud base is given by η = (z/zB )1/2 for all of the
updraft types where zB denotes the cloud base height.
49

3.2.4

Updraft Velocity and Entrainment Rate

The entrainment rate is deﬁned by

ϵ=

1 ∂Mu
Mu ∂z

(3.24)

and allowed to vary vertically. Based on the formulation of Gregory (2001), the updraft
velocity is calculated by
1 ∂ ŵ2= aB − ϵŵ2
2 ∂z

(3.25)

where w and B are the vertical velocity and the buoyancy of updraft air parcel respectively.
a is a dimensionless constant parameter ranging from 0 to 1 and represents a ratio of buoyancy
force used to accelerate the updraft velocity. The hats indicate the values of the updraft. The
second term on the right-hand side represents reduction in the upward momentum of the airparcel through the entrainment. Here and hereafter, the equation number corresponds to thatin Chikira and Sugiyama (2010).
Then it is assumed thatϵŵ2 ≃ Cϵ aB,(3.26)

where Cϵ is a dimensionless constant parameter ranging from 0 to 1. This formulation
denotes that a certain fraction of the buoyancy-generated energy is reduced by the entrainment,
which is identical to the fraction used to accelerate the entrained air to the updraft velocity.
Thus, the entrainment rate is written as

ϵ = CϵaB.
ŵ2(3.27)

Eqs. (3.25) and (3.27) lead to
1 ∂ ŵ2= a(1 − Cϵ )B
2 ∂z(3.28)

which shows that ŵ is continuously accelerated upward when buoyancy is positive. ManyCRM and LES results show, however, that updraft velocity is often reduced if the parcel approaches its cloud top. For this reason, adding an additional term, we use
1 ∂ ŵ2
1 ŵ2= a(1 − Cϵ )B −
2 ∂z
z0 2

(3.29)

where the last term denotes that the energy of the updraft velocity is relaxed to zero with aheight scale z0 . Eq. (3.29) is discretized as
1 ŵk+1/2 − ŵk−1/2
1 ŵk+1/2= a(1 − Cϵ )Bk −
2
∆zk
z0 2
2

2

2

50

(3.30)

where k is an index of full levels and k + 1/2 and k − 1/2 indicate the upper and lower sidesof the half levels. ∆z is the depth of the model layer. Note that the equation is solved for ŵ2
rather than ŵ.
The buoyancy of the cloud air parcel is determined by

B=
(
≃g

g
(T̂v − T̄v )T̄)ĥ − h̄∗+ ε(qˆv − q¯v ) − [(qˆl + qˆi ) − (q¯l + q¯i )]Cp T̄ (1 + γ)

(3.31)(3.32)

where g and Tv denote gravity and virtual temperature respectively. ε = Rv /Rd − 1 where
Rv and Rd are the gas constants for water vapor and dry air respectively.
ŵ, B and ϵ are calculated for each of the updraft types separately, but we omit the subscript
j for convenience.
3.2.5

Normalized Mass Flux and Updraft Properties

The properties of the updraft are determined by∂η ĥ= ϵη h̄ + Qi ,∂z

(3.33)∂η qˆt= ϵη q¯t − P∂z(3.34)

∂η
= ϵη,
∂z

(3.35)

andwhere Qi and P denote heating by liquid-ice transition and precipitation respectively. All
the other variables such as temperature, speciﬁc humidity, and liquid and ice cloud water arecomputed from these quantities. Tracers are calculated by a method identical to that for q̂t .
Equation (3.35) leads to∂ ln η
= ϵ.
∂z(3.36)ln ηk+1/2 − ln ηk−1/2= ϵk .
∆zk(3.37)

Then, η and ϵ are discretized as

Note that this discrete form leads to an exact solution if ϵ is vertically constant. Also, η is
ﬁnite as far as ϵ is. For ϵk , a maximum value of 4 × 10−3 m−1 is applied.
Equations (3.33) and (3.34) are written as51

∂η ĥ
= E h̄ + Qi ,∂z

(3.38)∂η q̂t= E q̄t − P∂z(3.39)

respectivuly, where E = ϵη. These equations are discretized asηk+1/2 ĥk+1/2 − ηk−1/2 ĥk−1/2= Ek h̄k + Qi,k∆zk

(3.40)

ηk+1/2 q̂t,k+1/2 − ηk−1/2 q̂t,k−1/2= Ek q̄t,k − Pk∆zk

(3.41)

Considering the relation that ∂η/∂z = ϵη, we descretize Ek as

Ek =ηk+1/2 − ηk−1/2∆zk(3.42)

Note that conservation of mass, energy, and water is guaranteed with Eqs. (3.37)–(3.42).
This set of equations leads to exact solutions of ĥ under the special case that ϵ and h̄ arevertically constant and Qi is zero. From Eqs. (3.37), (3.40), and (3.42), assuming Qi is zero,ĥk+1/2 = e−ϵk ∆zk ĥk−1/2 + (1 − e−ϵk ∆zk )h̄k ,(3.43)

which shows that ĥk+1/2 is a linear interpolation between ĥk−1/2 and h̄k . Thus, the stabilityof ĥ is guaranteed. The same property applies to q̂t as well, if P is zero.
These calculations are made for each of the updraft types separately, but we omit the subscript j for convenience.
3.2.6

Spectral Representation

Following the spirit of the Arakawa–Schubert scheme, updraft types are spectrally represented.
Diﬀerent values of cloud-base updraft velocities are given from the minimum to the maximum
values with a ﬁxed interval. The minimum and maximum values are set to 0.1 and 1.4 ms−1 ,with an interval of 0.1 ms−1 for MIROC6. The minimum and maximum values and the interval
can be changed.
Then, the updraft properties are calculated upward with Eqs. (3.27), (3.29), (3.33), (3.34),
and (3.35). This upward calculation continues even if the buoyancy is negative as long as
the updraft velocity is positive. If the velocity becomes negative at some level, the air parcel
detrains at the neutral buoyancy level which is below and closest to the level. That is, the scheme
automatically judges whether the rising parcel can penetrate the negative buoyancy layers when
there is a positive buoyancy layer above. The eﬀect of the convective inhibition (CIN) near cloudbase is also represented by this method. Note, however, that an eﬀect of overshooting above
cloud top is not represented for simplicity (i.e., detrainment never occurs above cloud top).

52

3.2.7

Cloud-Base Mass Flux

The cloud-base mass ﬂux is determined with the prognostic convective kinetic energy closure
proposed by Arakawa and Xu (1990). That is, the cloud kinetic energy for each of the updraft
types is explicitly predicted byK∂K= AMB −,∂tτp(3.44)

where K and A are the cloud kinetic energy and cloud work function respectively, and τpdenotes a time scale of dissipation. The cloud work function A is deﬁned asA≡

Z zTBη dz .

(3.45)zB

The cloud kinetic energy is linked with MB by
K = αMB2 .

(3.46)

The cloud-base mass ﬂux is then solved for each of the updraft types.
3.2.8

Microphysics

The method to obtain temperature and speciﬁc humidity of in-cloud air from moist static energy
is identical to that in Arakawa and Schubert (1974). The ratio of precipitation to the total
amount of condensates generated from cloud base to a given height z is formulated as
Fp (z) = 1 − e−(z−zB −z0 )/zp ,(3.47)

where z0 and zp are tuning parameters.
The ratio of cloud ice to cloud condensate is determined simply by a linear function oftemperature,1
Fi (T ) = (T2 − T )/(T2 − T1 )0T ≤ T1T1 < T < T2
T ≥ T2(3.48)

where T1 and T2 are set to 258.15 and 273.15 K. The ratio of snowfall to precipitation is also
determined by this function.
From the conservation of condensate static energy, Cp T + gz + Lv q − Li qi where Li is thelatent heat of fusion, for a cloud parcel, Qi in Eq. (3.33) is written asQi = Li∂η q̂i− ϵη q̄i∂z(3.49)and discretized asQi k = Liηk+1/2 q̂i,k+1/2 − ηk−1/2 q̂i,k−1/2− Ek q̄i,k
∆zk
53


(3.50)

Strictly, the ratio of the cloud ice to the cloud condensate should be recalculated after themodiﬁcation of temperature by Qi and the iterations of the calculation are required; however,
it is omitted for simplicity.
Melting and freezing of precipitation occurs depending on wet-bulb temperature of large-scale
environment and cumulus mass ﬂux.
3.2.9

Evaporation, Sublimation and Downdraft

A part of precipitation is evaporated at each level asEv = ae (q̄w − q̄)

P
VT


,

(3.51)

where Ev , qw and VT are the mass of evaporation per a unit volume and time, wet-bulb
saturated speciﬁc humidity and terminal velocity of precipitation respectively ae is a constant.
Downdraft mass ﬂux Md is generated as∂Md= −be ρ̄(T̄w − T̄ )P,∂z(3.52)

where ρ and Tw are density and wet-bulb temperature, respectively; be is a constant. Properties of downdraft air are determined by budget equations and the detrainment occurs at neutral
buoyancy level and below cloud base.
If the precipitation is composed of both rain and snow, the rain (snow) is evaporated (sublimated) in the same ratio as the ratio of rain (snow) to the total precipitation when the precipitation evaporates to produce downdrafts.
3.2.10

Cloudiness

Fractional cloudiness of the updrafts Cu used in the radiation scheme is diagnosed byCu =XCmax − Cmin(lnMu,j − ln Mmin ) + Cmin ,ln Mmax − ln Mmin

(3.53)jwhere Cmax , Cmin , Mmax , Mmin are the maximum and minimum values of the cloudiness andcumulus mass ﬂux respectively.
The grid mean liquid cloud mixing ratio in the updrafts is given by

lc =

βCu X
q̂l,j Mu,j ,M(3.54)jwhere β is a dimensionless constant. The grid mean ice cloud mixing ratio is determined
similarly.
3.2.11

Cumulus Momentum Transport

Following Gregory et al. (1997), the zonal and meridional velocities of the updrafts are calculatedas54∂η û∂ ū
= ϵηū + Cm η ,∂z∂z(3.55)

where Cm is a constant from 0 to 1 representing the eﬀect of pressure. This equation can berewitten as
∂η û
∂η ū
= (1 − Cm )ϵηū + Cm,∂z∂z

(3.56)

ηk+1/2 ūk+1/2 − ηk−1/2 ūk−1/2
ηk+1/2 ûk+1/2 − ηk−1/2 ûk−1/2= (1 − Cm )Ek ūk + Cm
.
∆zk
∆zk

(3.57)and is discretized as

The horizontal velocities of the downdrafts are calculated similarly. The tendencies of zonal
and meridional velocities by the cumulus momentum transport (CMT) are calculated as


∂u
∂t





∂v
∂t

= −g
CMT,k

(ρu′ w′ )k+1/2 − (ρu′ w′ )k−1/2,∆pk(3.58)(ρv ′ w′ )k+1/2 − (ρv ′ w′ )k−1/2∆pk(3.59)


= −g
CMT,k

respectively, where ρu′ w′ and ρv ′ w′ are total zonal and meridional momentum ﬂuxes respectively and ∆pk = pk − pk+1 .

55

3.3
3.3.1

Shallow Convection Scheme
Overview of Shallow Convection

Shallow convection is the most frequent type of convective cloud in the tropics and subtropics
Its impact on climate through the energy budget due to atmospheric radiation is considered
important (Stevens, 2005). Shallow convection is responsible for transporting the boundary layer
air to the free atmosphere. It is often not accompanied by precipitation and is characterized by
the fact that precipitation-induced downdraft does not reach the surface as in deep convection.
This section brieﬂy describes the vertical structure of the boundary layer favorable for shallow
convection. When the ground surface is heated by sunlight or cold air ﬂows in from above, theenergy of convective instability is dissipated by turbulence in the bottom of the atmosphere,
forming a mixed layer with a nearly uniform vertical distribution of potential temperature andwater vapor at a thickness of about 600 m to 800 m from the surface. At the upper end of
the mixed layer, there is a transition layer of weakly stable stratiﬁcation, which is the heightat which water vapor in updraft begins to condense (lifting condensation level, LCL). Above
LCL, the temperature decreases according to the moist adiabatic lapse rate, and the updraft
is observed as clouds. Above the level of free convection (LFC), the cloud continues to growwhile mixing with surrounding air. The growth of these convective clouds is limited by thetemperature inversion layer at the lower end of the free atmosphere, and the cloud tops areoften located about 2 km from the surface.
In the former versions of MIROC, A cumulus parameterization proposed by Chikira andSugiyama (2010) deals with multiple cloud types including shallow cumulus and deep convective clouds. However, it tends to overestimate low-level cloud amounts. To cope withthis bias and improve the performance for reproducing current climate, the shallow convection scheme is introduced from the 6th version of MIROC (Tatebe et al., 2019, Ogura et al.,2017, Ogura, 2015)．The source code in concern (pshcn.F) consists of SUBROUTINE:[PSHCN]
and SUBROUTINE:[DISTANCE]. The input values for SUBROUTINE:[PSHCN] are temperature, water vapor mixing ratio, and liquid water mixing ratio, ice mixing ratio. It predicts liquid water
potential temperature, total water mixing ratio, ice mixing ratio, and horizontal componentsof wind in response to vertical transport. It also diagnoses cloud fraction and precipitation.
SUBROUTINE:[DISTANCE], which is called inside SUBROUTINE:[PSHCN], calculates the degree of
buoyancy-induced updraft and mixing with the environment. Since the variables diagnosedin the cumulus scheme (SUBROUTINE:[CUMLUS]) are referenced to determine the conditions forshallow convection, SUBROUTINE:[PSHCN] is required to be run after SUBROUTINE:[CUMULUS],
followed by the diagnosis of cloud fraction. On the other hand, it should be run before the land
surface process SUBROUTINE:[SURFCE] because precipitation by convection is referenced in theland surface and ocean models.
3.3.2

Basics of Cloud Model

Subgrid clouds are modeled based on the frameworks proposed by Bretherton et al. (2004) andPark and Bretherton (2009). This scheme employs a simple plume model for cloud to calculatevertical transport of conserved variables and precipitation due to updraft. An ensemble ofshallow convection in a horizontal grid, which is expressed as a single updraft plume, is supposed
to experience horizontal mixing with the environment (entrainment/detrainment). The ﬂux ofvertical transport of mass is assumed in the following form:
ρw′ ψ ′ ≈ Mu (ψu − ψ),(3.60)

where Mu = ρu σu wu is mass ﬂux of updraft (ρu ,σu , and wu stand for density in updraft, area
fraction of updraft in a grid, and vertical velocity, respectively)，ψu is a conserved variable56

transported by convection (e.g. liquid water potential temperature, total water mixing ratio,
horizontal components of momentum) in updraft，ψ denotes the average value in the environmental ﬁeld of the same conserved value. The eﬀects of vertical transport due to shallow
convection are represented by determining the vertical proﬁles of unknown values Mu and ψu .
Flux of mass and conserved values are diagnosed as
∂Mu
=E−D
∂z

(3.61)∂(ψu Mu ) = Xψ + Sψ Mu ,∂z

(3.62)

where Xψ represents horizontal mixing with environmental air, and Sψ is source term. E and
D are rates of entrainment and detrainment, which are described in fractional fromE = ϵMu(3.63)

D = δMu .

(3.64)

Substituting ψ for grid value and assuming the horizontal mixing term as Xψ = Eψ − Dψuresults in∂Mu= Mu (ϵ − δ)∂z∂ψu
= ϵ(ψ − ψu ) + Sψ .
∂z

(3.65)(3.66)In MIROC6, changes in liquid water potential temperature due to precipitation and the eﬀect ofsubgrid pressure gradient on horizontal momentum are included in Sψ . Consequently, equations(3.65) and (3.66) results in a closure problem of two parameters δ and ϵ. By determining δ andϵ by the formulation described in section 3.3.3.4 and solving diﬀerential equations along withboundary condition at cloud base, vertical proﬁles of Mu and ψu are calculated.
3.3.3

Computation in PSHCN

The eﬀect of convective updraft is calculated as follows.
• Liquid water potential temperature θl and total water qt are diagnosed from input temperature T , water vapor mixing ratio qv , liquid water mixing ratio ql , ice mixing ratioqi ,• Updraft mass ﬂux at cloud base is diagnosed.
• Height of cloud base is diagnosed.
• Presence of shallow convection is determined.
• Vertical proﬁles of Mu , θl , qt , horizontal wind components u and v are diagnosed.
• θl , qt , qi , u, v, liquid water temperature Tl are predicted.
• T , qv , and ql are diagnosed according to Tl and qt .

57

3.3.3.1

Lower boundary condition: diagnosis of cloud base mass ﬂux

The mass ﬂux at cloud base is formulated as it depends on turbulent kinetic energy (TKE)in boundary layer and convective inhibition (CIN) at the top of boundary layer.
Firstly, the vertical proﬁle of updraft velocity is supposed to fulﬁll
1 ∂ 2
w = aBu − bϵwu2
2 ∂z u(3.67)

all over the layers with shallow convection. Bu means updraft buoyancy, a and b are empiricalparameters. The ﬁrst term of the right-hand side of (3.67) is acceleration by buoyancy, andthe second term represents drag by entrainment. By assuming no entrainment below LFC andintegrating (3.67) from cloud base to LFC, The critical value of upward velocity for updraft
plume to reach LFC, wc , can be determinedwc =p2a(CIN ).

(3.68)

Updrafts that exceed this critical value penetrates from cloud base.
Computation of CIN is based on Appendix C of Bretherton et al.,
CIN = [Bu (pbase ) + Bu (pLCL )]pLF C − pLCL
pLCL − pbase+ Bu (pLCL ).
g(ρLCL + ρbase )g(ρLF C + ρLCL )(3.69)In the following, subscript base represents the value at the top of mixing layer. In MIROC6, forsimplicity, Bu (pbase ) is set to zero.
Secondly, to obtain the information of vertical velocity at cloud base, the statistical distribution of w is assumed to follow Gaussian distribution

1w2
f (w) =
exp −(3.70)
2πkf eavg
2kf eavg
with variance equal to kf eavg , where eavg is average TKE diagnosed in turbulent and vertical
diﬀusion scheme. kf is an empirical parameter describing the partitioning of TKE between
horizontal and vertical motions at the subcloud layer inversion, whose recommended value basedon large eddy simulation is 0.5.
By taking average of vertical velocity above the critical value wc , cloud base mass ﬂux Mu,baseis diagnosed as
r


Z ∞
kf eavg
wc2
Mu,base = ρbase
wf (w)dw = ρbaseexp −,(3.71)
2π
2kf eavg
wc
where ρbase is density at LFC. This mass ﬂux is larger for larger boundary layer TKE and smaller
for larger CIN.
3.3.3.2

Diagnosing height of cloud base

The cloud base height is set between the top of the boundary layer and the LCL. The larger
the CIN is, the lower the cloud base becomes. The top of boundary layer is diagnosed as thelevel with maximum vertical gradient of relative humidity. Let zHi be the higher of this leveland LCL, and zLo be the lower, then the cloud base altitude zbase is set
zbase = zHi − (zHi − zLo )CIN − CINLo.
CINHi − CINLo(3.72)

CINHi and CINLo are coeﬃcients which satisfy CINLo ≤ CIN ≤ CINHi for a typical value ofCIN.
58

3.3.3.3

Determination of the presence of shallow convection

For each horizontal column, whether shallow convection occurs is determined with following
criteria.
• If estimated inversion strength (EIS; Wood and Bretherton, 2006) exceeds a certain threshold, the environmental ﬁeld is judged to be dominated by stratocumulus clouds, and shallow convection is not generated. This criterion is introduced because the vertical resolutionof climate models does not suﬃciently represent the thin and strong inversion layer overthe boundary layer, and underestimates CIN, which leads to an overestimation of shallow
convection. EIS is estimated by EIS = θ700 − θ0 − Γ850
m (z700 − LCL) where θ700 and θ0 are
potential temperature at 700hPa and surface, Γ850
ismoist adiabatic lapse rate at 850hPa,m
and z700 is height of 700hPa.
• If the intensity of cumulus convection diagnosed by SUBROUTINE:[CUMULUS] exceeds athreshold, the environmental ﬁeld is supposed to be dominated by deep convection andshallow convection is not generated.
• If the areal fraction of shallow convection is under a threshold, computation of shallow
convection is omitted.
3.3.3.4

Diagnosing vertical proﬁle of updraft mass ﬂux

For the grid boxes that contain shallow convection, entrainment and detrainment is calculated using the value of ψu at cloud base and Mu,base . Fractional entrainment and detrainment
are computed based on the framework of buoyancy sorting suggested by Kain and Fritsch (1990).
In a layer of thickness δz, equal parts ϵ0 Mu δz of updraft and environmental air are involved in
the lateral mixing process that creates a spectrum of mixtures. This yields a total mixing mass
ﬂux 2ϵ0 Mu δz, with fractional mixing rate ϵ0 = c0 /H (c0 is a certain empirical coeﬃcient and H isthe height from surface). In the mixed air, there exists states with probability density such thatthe air from the environmental ﬁeld occupies a proportion χ. Here, for simplicity of calculation,
it is considered that the state from pure moist air (χ = 0) to pure environmental air (χ = 1)
is distributed with uniform probability (Kain-Fritsch scheme assumes Gaussian distribution).
Based on the buoyancy force on the mixed air, the entrainment or detrainment is determined.
SUBROUTINE:[DISTANCE] is called in SUBROUTINE:[PSHCN]. The output variables in this subroutine are liquid water potential temperature (THETLU) and bool value for entrainment or
detrainment (JUDGE).
The occurrence of entrainment is judged as follows. Firstly, if the updraft air is not saturated,
entrainment is not assumed to occur. Nextly, with virtual potential energy in the environmentalﬁeld (θv ) and updraft (θvu ), buoyancy force on the parcel is deﬁned:
Bu = gθvu − θv
θv

(3.73)and entrainment occurs when the buoyancy on parcel is positive. Furthermore, even when thebuoyancy is negative, entrainment occurs if the parcel can travel longer than a certain eddy
mixing distance lc = c1 H, where c1 = 0.1 is an empirical constant, chosen to optimize thetrade-cumulus case. This criterion corresponds to the critical buoyancy value
Bc = −

1 wu2
2 lc(3.74)and otherwise, all the mixed air is detrained. Therefore, Once the critical value of the mixing
state χc is obtained, which allows the updraft to rise a distance lc under negative buoyancy,
59

the air in the environmental ﬁeld entrained into the cloud and the air in the updraft that isdetrained can be determined as follows
Z χc
Mu ϵ = 2ϵ0 Mu
χq(χ)dχ = ϵ0 Mu χ2c(3.75)0
Z 1
Mu δ = 2ϵ0 Mu(1 − χ)q(χ)dχ = ϵ0 Mu (1 − χc )2 .
(3.76)
χcThus, letting
ϵ = ϵ0 χ2c(3.77)

δ = ϵ0 (1 − χc )2 ,(3.78)

equatinons (3.65) and (3.66) are expressed as follows
1 ∂Mu= ϵ − δ = ϵ0 (2χc − 1)
Mu ∂z
∂ψu
= ϵ(ψ − ψu ) + Sψ = ϵ0 χ2c (ψ − ψu ) + Sψ ,∂zwhere χc is computed based on virtual potential temperature of mixed air




βL
θv (χ) = θvu + χ β(θl − θlu ) −− θu (qt − qtu )cp Π(3.79)(3.80)(3.81)(Bretherton et al., 2004). β is a thermodynamic parameter which depends on temperature andpressure deﬁned by Randall (1980), θlu is liquid water potential temperature in updraft, θu is
updraft potential temperature，qt is total water mixing ratio of environment，qtu is total water
mixing ratio of updraft, L is latent heat of vaporization，cp is speciﬁc heat capacity of dry air
at constant pressure, and Π is the Exner function.
Consequently, the governing equations (3.67), (3.79), and (3.80) for vertical proﬁles of wu ,Mu , and ψu are obtained. These equations are discretized and integrated upward one layer at
a time using the lower boundary condition in section 3.3.3.1 to yield the vertical proﬁle of eachvariables.
Afterward, from liquid water potential temperature and total water mixing ratio, liquid water
mixing ratio ql and water vapor mixing ratio qv are diagnosed. The cloud water that exceeds
a threshold is disposed as rainwater qr , and liquid water potential temperature is updated
according to the amount of qr . This corresponds to Sψ in (3.81).
The formulation of the vertical ﬂux in this scheme is equal to the assumption that the updraft
is not large enough to replace all of the air in a grid box in the time step ∆t. Therefore, thefollowing limiter is imposed to prevent numerical instability when diagnosing mass ﬂux of the
updraft.
ρ∆zMu = min. Mu ,(3.82)
∆t

60

3.4

Large Scale Condensation

The SUBROUTINE:[PDF2CLD] and SUBROUTINE:[CLD2PDF] are written in pmlsc.F ﬁle. These
are called in padmn.F, pcumc.F, pshcn.F, pcldphys.F and pvdfm.F ﬁles.
3.4.1

Physical Basis for Statistical PDF Scheme

General Circulation Models (GCMs) typically adopt fractional cloud cover (the volume of cloudy
air per total air volume in a grid box) assumption to realistically represent clouds because of
their coarse horizontal resolution (O(100km)). Statistical cloud schemes assume a subgrid‐scale
probability density function (PDF) of humidity within the grid. Integration of the speciﬁc PDFs
will give the cloud fraction and the amount of water condensate consistently.
By means of the “fast condensation” assumption, the cloud water amount in a local area in
the grid is
qc = (qt − qs ) δ (qt − qs )(3.83)

where qs denotes the saturation //mixing ratio and qc does the cloud water mixing ratio. qt
is sum of water vapor and cloud water mixing ratio. δ(x) denotes the Heviside function of x.
The majority of statistical cloud schemes use the so-called “s-distribution” following Sommeria and Deardorﬀ (1977). A single variable s, which considers the subgrid-scale perturbationsof liquid temperature Tl and total water mixing ratio qt , is employed. s is deﬁned ass = aL (qt − αL Tl )(3.84)

aL = 1/ (1 + LαL /cp ) , αL = ∂qs / ∂T |T =T̄l .

(3.85)

where

For any choice of the PDF of s, denoted as G(s), the grid-mean cloud fraction, C, and cloudwater content, qc , are obtained by integrating G(s) and (Qc + s)G(s),Z ∞G(s)ds(3.86)(Qc + s) G(s)ds,(3.87)

C=
−Qc

Z ∞
q̄c =

−Qc

where Qc denotes the grid-scale saturation deﬁcit deﬁned as
Qc ≡ aL q̄t − qs T̄l , p̄ .
3.4.2

(3.88)

Hybrid Prognostic Cloud scheme

The statistical scheme implemented in MIROC6 is called Hybrid Prognostic Cloud (HPC)
scheme (Watanabe et al. 2009). The HPC scheme proposes two types of shape for the PDF
G(s), Double-uniform PDF and Skewed-triangular PDF. Here we focus on Skewed-triangular
scheme because MIROC6 adoptes the shape. The physical basics of the scheme are in common
with Double-uniform PDF.
61

Example of the basis PDF for HPC: skewed-triangular functions. Copied from Fig.1 in
Watanabe et al. 2009.
The scheme predicts variance (V ) and skewness (S) of the PDF. V , S, the second momentµ2 , and the third moment µ3 are deﬁned as follows.
µ2 ≡ V =3/2

Z ∞s2 G(s)ds(3.89)−∞µ3 ≡ µ2 S =Z ∞s3 G(s)ds(3.90)−∞V and S are aﬀected by cumulus convection, cloud microphysics, turbulent mixing, andadvection.
The integrals to obtain C and qc are symbolically expressed as
C = IC p̄, T̄l , q̄t , V, Sq̄c = Iq p̄, T̄l , q̄t , V, S



(3.91)(3.92)

where p̄ denotes the pressure. The overbars denote the grid-mean quantity.
If the PDF is not too complicated, (1, 2) can be analytically solved for V and S by deﬁning
integrand functions I˜ as
V = I˜V p̄, T̄l , q̄v , q̄c , C
S = I˜S p̄, T̄l , q̄v , q̄c , C(3.93)(3.94)The relationship between (1, 2) and (4, 5) is quasireversible. The double-uniform function and skewed-triangular function PDFs are selected for G(s) because of their feasibility in˜
analytically nalderiving I.
3.4.3

PDF Change Through Processes

The HPC cloud scheme is composed using prognostic equations for four variables determining
˜I, namely, Tl , qt , V , and S. The prognostic variables can be Tl , qt , C, and qc that determine I.
Prognostic equations for the PDF variance and skewness are expressed as
DV
∆V
∆V
∆V
∆V
=
+
+
+
− εV
Dt
∆t conv.
∆t micro.
∆t turb.
∆t others(3.95)

∆S
∆S
∆S
DS
∆S
+
+
+
− εS
=
Dt
∆t conv.
∆t micro.
∆t turb.
∆t others(3.96)

where subscripts ‘conv.’, ‘micro.’ and ‘turb.’ indicate cumulus convection, cloud microphysics and turbulent mixing processes respectively, which all aﬀect the PDF shape. The lastterms represent dissipation due to subgrid-scale horizontal motions. The speciﬁc formulations
for each term are described below.
The HPC scheme is referred to as and G(s) is updated every after the process that aﬀects
cloud water PDF. G(s) is thus modiﬁed several times within a single time step.
62

3.4.3.1

Cumulus Convection

The total eﬀect of cumulus convection to the PDF moments is written as
∆V
∂V
∆I˜V
= Mc
+
∆t conv.
∂z∆t(3.97)

∆S
∂S
∆I˜S
= Mc
+
∆t conv.
∂z∆t(3.98)

Mc is the cumulus mass-ﬂux including updraft in the convection tower and downdraft in theenvironment. The vertical transport of the PDF moments is represented by the ﬁrst terms on
the right side hand of (14, 15).
Cumulus convections modify the grid-mean Tl , qt , and qc by upward transportation of gridmean moist static energy, qv , and qc . Detrainment also aﬀects these variables. The detrainmentof the cloudy air mass is considered, as in Bushell et al. (2003),∂C= D(1 − C)∂t conv.

(3.99)The second terms on the right hand side of (14, 15) indicates that the changes in the PDF
moments is calculated consistent with the changes in the grid-scale temperature, humidity, cloudwater, and cloud fraction.
∆I˜X = I˜X p̄, T̄l + ∆T̄l , q̄v + ∆q̄v , q̄c + ∆q̄c , C + ∆C − I˜X p̄, T̄l, q̄v , q̄c , C(3.100)(3.101)
where X is either V or S.
3.4.3.2

Cloud Microphysics

The tendency due to microphysical processes can be written in a similar manner to the
cumulus convection eﬀect.
∆V
∆I˜V
=
∆t micro.
∆t(3.102)∆S∆I˜S=∆t micro.
∆t(3.103)

Changes in T̄l , q̄v , andq̄c are derived from microphysical tendency terms including precipitation, evaporation,and melting/freezing.

63

3.4.3.3

Turbulent Mixing

From the deﬁnition of s, the PDF variance V becomes2
Πθl′2 − 2αL Πqt′ θl′ ,V = a2L qt′2 + αL

(3.104)

where Π is the Exner function. Assuming the level-2 closure in Nakanishi and Niino (2004),
the time evolution of V can be derived as
"

2


∂ θ̄l
∂ q̄t 2∆V
2
=2aL (αL Π)2 KH+ Kq∆t turb.
∂z∂z∂ θ̄l ∂ q̄t
2q−αL Π (KH + Kq )
−
V,
∂z ∂z
Λ2

(3.105)where KH and Kq are the mixing coeﬃcients for sensible heat and moisture, respectively.
2
q = u′2 + v ′2 + w′2 denotes the turbulent kinetic energy. The other symbols follow the original
notation.
Since the turbulence production does not aﬀect the PDF shape parameter deﬁned by thethird moment (cf. Tompkins 2002), the skewness change ∆S/ ∆t|turb. is simply calculated dueto the variance change in (28).
3.4.3.4

Subgrid-Scale Horizontal Eddy

In the planetary boundary layer, the subgrid-scale inhomogeneity is dissipated due to theturbulent mixing. In free atmosphere, the grid box will be homogenized mainly due to mesoscale
motions, which are expressed by the Newtonian damping as in (Tompkins 2002): εV = τVh , εS =Sτh , where the relaxation timescale is parameterized by the horizontal wind shear asτh−1 = Cs2

(

∂ ū
∂x

2


+

∂v̄
∂y

2 )1/2(3.106)The coeﬃcient Cs is set to 0.23 following Tompkins (2002).
3.4.3.5

Other Processes

Dynamics, shallow convection, radiation, mass source, and dissipation heating processeschange the grid-mean temperature and humidity. Such eﬀects on the shape of PDF are included
following (16).
3.4.4Solving Procedures

The shape of the Skewed-triangular PDF is represented as follows. The widths deﬁned bypositions of the left and right edges on the s-coordinate are denoted as a and b, respectively.
The position of the top, denoted as q, is constrained by a + b + q = 0. By deﬁnition, q ≤ b anda ≤ q must be satisﬁed. The PDF is then expressed as(
G(s) =

2(s−b)− (b−q)(b−a)

for q < s ≤ b2(s−a)(q−a)(b−a)

for a < s ≤ q

64(3.107)

The pmlsc module includes two main subroutines, PDF2CLD and CLD2PDF. The subroutine PDF2CLD calculates C and q̄c given p̄, Tl, q̄t , V, S. The subroutine CLD2PDF calculates V
and S given p̄, Tl, q̄t , q̄c , C. We will derive the concrete calculation processes in this subsection.
3.4.4.1

Caluculation of Cloud Variables from PDF Moments

This is written in SUBROUTINE:[PDF2CLD].
From µ1 , µ2 , µ3 to a, b, q
The ﬁrst, second, and third moments of the PDF is calculated as follows.
Z q+b
sG(s)ds

µ1 =

b−a
3

(3.108)a2 + ab + b2
18

(3.109)

=q+

q−a

Z q+b
µ2 =

(s − µ1 )2 G(s)ds =

q−a

Z q+b

(b − a) 2a2 + 5ab + 2b2µ3 =(s − µ1 ) G(s)ds =
270
q−a



3

(3.110)From (7,8,9), we will derive the solution for a, b, q given µ1 , µ2 , µ3 .
We deﬁne δ ≡ b − a, β ≡ ab. (8,9) areδ 2 + 3β = 18µ2

(3.111)δ (β + 12µ2 ) = 90µ3

(3.112)Eliminate β or δ from these equations, you will get the equations.
δ 3 − 54µ2 δ + 270µ3 = 0(3.113)1
β = 6µ2 − δ 2
3(3.114)

We apply the formula for the solution of a cubic equation to (10) to obtain δ.




p
1
−135µ3  4 
δ = 2 18µ2 cos  cos−1  q+ π
3
3(18µ2 )3
β is obtained from (11). We deﬁne α ≡as follows.

(3.115)pδ 2 + 4β for simplicity. Finally, a, b, q is calculated65a = (α − δ)/2

(3.116)b = (α + δ)/2

(3.117)q = µ1 − δ/3(3.118)From PDF to C and qc
Once the PDF G(s) is determined by the parameters a, b, q, the cloud fraction C and gridmean cloud water mixing ratio q¯c are derived as follows.

C=



0



 (Qc +b)2

if q ≤ −Qc ≤ b





 1

if a ≤ −Qc ≤ q
if − Qc < a(b−q)(b−a)(Qc +a)2(q−a)(b−a)

if b < −Qc



 01

q¯c =

3.4.4.2

3 C (Qc + b)
Q
− 1 (1 − C) (Qc + a)


 c 3
Qc

if b < −Qc
if q ≤ −Qc ≤ b
if a ≤ −Qc ≤ q
if − Qc < a(3.119)(3.120)Caluculation of PDF moments from cloud variables.

This is written in SUBROUTINE:[CLD2PDF]
From q¯c , C to a, b, q
We can not determine the position of Qc in the triangle at the beginning of the calculation.
Thus we calculate a, b assuming that a ≤ −Qc ≤ q at ﬁrst. If the calculated parameters are
physically consistent with the PDF (a + b ≥ 0), a, b, q are determined. Otherwise, we regardq ≤ −Qc ≤ b and then a, b, q are derived.
1. a ≤ −Qc ≤ q
From (16), a is derived as follows.

a=3 (Qc − qc )
− Qc
1−C

(3.121)

We eliminate q from (15) using q = −a − b. The quadratic equation for b is obtained.
b2 + ab − 2a2 + (Qc + a)2 /(1 − C) = 0(3.122)The physically meaningful solution for b is

b=

 q

−a 9a2 − 4 (Qc + a)2 /(1 − C) /2
66(3.123)

2. q ≤ −Qc ≤ b
From (16), b is
b=

3qc
− Qc
C

(3.124)

We eliminate q from (15) using q = −a − b. The quadratic equation of a is obtained.
a2 + ab − 2b2 + (Qc + b)2 /C = 0(3.125)The physically meaningful solution for a is


q
2
2
a = −b − 9b − 4 (Qc + b) /C /2(3.126)

Adjustment of Cloud Fraction
When there is no physically meaningful solution for (18), C is adjusted so that a reasonable
solution is obtained. The critical conditions for the existence of real solutions for (18) are as
follows.
9a2 − 4 (Qc + a)2 /(1 − C) = 0 (a ≤ −Qc ≤ q)
9b2 − 4 (Qc + b)2 /C = 0(q ≤ −Qc ≤ b)

(3.127)

Eliminate a and b uging (17), we get the relationship between C and qc ,


2

2
3(Qc −qc )3(Qc −qc )
4
−
Q
=
c
1−C
1−C
1−C

2
 2
3qc
3q
9 C − Qc = C4 Cc

9(a ≤ −Qc ≤ q)We
√ take the square root of the both sides of the equations and deﬁne γ1 ≡γ2 ≡ C. The cubic equations for γ is obtained.

γ13 − 3 1 − Qqcc γ1 ± 2 1 − Qqcc = 0 (a ≤ −Qc ≤ q)
γ23 − 3 Qqcc γ2 ± 2 Qqcc = 0(3.128)(q ≤ −Qc ≤ b)(q ≤ −Qc ≤ b)√

1 − C and(3.129)We deﬁne R1 = 1 − Qqcc , R2 = Qqcc .
 −4R sinh2 1 sinh−1 √ 1(R < 0) 3
  −R 
γ2 = 4R cos2 1 cos−1 √1 + 4 π(R > 1)
3
3
R(3.130)Note that, γ = γ1 , R = R1 (a ≤ −Qc ≤ q) or γ = γ2 , R = R2 (q ≤ −Qc ≤ b).
The actual calculation procedure is as follows. If the solution for (18) is not a real number,
C is adjusted using (26). Then we solve (18) again.

67

From a, b, q to µ2 , µ3
By deﬁnition, the PDF moments are expressed in terms of a and b.

3.4.4.3µ2 =a2 + ab + b2
6

(3.131)µ3 =−(a + b)ab10

(3.132)Treatment of Cloud Ice and in-Cloud Water Vapor

Because the original HPC scheme by Watanabe et al. (2009) does not consider the cloud
ice, it is modiﬁed when coupled with the Wilson and Ballard (1999) ice microphysics. Since thestatistical PDF scheme employs a ‘fast condensation’ assumption that is no more valid for ice,
the ice mixing ratio is assumed to be conserved in the large scale condensation process.
Here we assume that - the water vapor mixing ratio within the cloudy area in a grid isconstant - cloud ice preferentially exists in areas with large total water content
Based on these assumptions, the cloud fraction and each condensate mixing ratios are diagnosed. The notations for the mixing ratios (ql , qi , qv , qvi ) of liquid water (subscript l), ice(subscript i), vapor (subscript v), in-cloud vapor (subscript vi) are employed.
At ﬁrst the total condensate mixing ratio qc = ql + qi is diagnosed from qt and Tl assumingthat ice does not exist in the grid. The saturation mixing ratio is set for liquid (qsatl ).
Mixed-phase cloud is generated when the condensate amount is more than the ice content(qc > qi ), whereas the cloud fraction and vapor amount are adjusted in the case of a pure ice
cloud when the condensate amount is less than the ice content (qc < qi ). Speciﬁcally, qc , C andqvi are calculated as follows.
1. qc > qi
Liquid-phase clouds and ice clouds coexist.
ql = qc − qi

(3.133)

qvi = qsatl(3.134)

2. qc < qi
Only ice clouds exist (ql = 0). In this case, C and qvi are rediagnosed. We eliminate Qc in(15,16) assuming that qc = qi . Equations for C are given as

C3 =

9qi2
(q ≤ −Qc ≤ b)(b − q)(b − q)

C 3 + 3C 2 = 4 −9 (qi + a)2(a ≤ −Qc < q)
(q − a)(b − a)68

(3.135)(3.136)

From these equations, C is obtained as follows.
 q


9qi2
(b−q)2
3

0
≤
q
≤

i

(b−q)(b−a)
3(b−a) 9(qi +a)2(b−q)2
1−1C=
1 − 2(q−a)(b−a)−12 cos 3 cos
3(b−a) < qi ≤ −a1(−a < qi )(3.137),whereQc =p3qi− b = 3 3qi (b − q)(b − a) − b.
C

(3.138)Given Qc , qvi = qt − Qc is calculated as follows.

3qi


 qt − C + b
qvi =



2

(b−q)
0 ≤ qi ≤ 3(b−a)

3(qi +a)(b−q)2
q
−
+
a
<
q
≤
−a
t
i
2+C
3(b−a)qt − qi(−a < qi )69

(3.139)

3.5

Cloud Microphysics

The SUBROUTINE:[CLDPHYS] is written in the pcldphys.F ﬁle.
3.5.1

Overview of Cloud Microphysics

Cloud microphysics control the conversion from water condensate to precipitate. The condensate
parameterization closely links to the lifetime of and radiative properties of the clouds.
The stratiform (non-convective) cloud microphysics in MIROC6 (Tatebe et al. 2019) arebasically the same as those used in MIROC5 (Watanabe et al. 2010). MIROC5 implemented aphysically based bulk microphysical scheme. The previous version of the scheme in MIROC3.2
diagnoses the fraction of liquid-phase condensate to total condensate simply as a function of
the local temperature. In contrast, the explicit treatment of ice cloud processes allows ﬂexible
representation of the cloud liquid/ice partitioning in MIROC5 and MIROC6 (Watanabe etal. 2010; Cesana et al. 2015).
The MIROC6 cloud microphysics scheme uses four quantities to describe water in the atmosphere: vapour; liquid-phase cloud droplets; raindrops; and frozen water. Only one quantity,
which we will refer to as ‘ice’, is used to describe all frozen water in large-scale clouds, including
aggregated snow, pristine ice crystals and rimed particles. Physically based transfer terms link
the four water quantities. The scheme treats two prognostic condansate variables: ice water
mixing ratio qi and cloud water mixing ratio qc . Water vapor mixing ratio qv aﬀects the rate of
microphysical processes and qv itself is also modiﬁed via microphysical processes. Ice number
concentration Ni is diagnosed as a function of qi and air temperature T in K. Cloud number
concentration Nc is predicted by the online aerosol module implemented. Rain water mixing
ratio qr is treated as a diagnostic variable: qr falls out to the surface within the time step. Cloud
fraction is predicted as described in the section ‘pmlsc: Large Scale Condensation’.
The cold rain parameterization following Wilson and Ballard (1999) predicts qi using physically based tendency terms, which represent homogeneous nucleation, heterogeneous nucleation,
deposition/sublimation between vapor and ice, riming (cloud liquid water collection by fallingice), and ice melting. The warm rain processes produce rain as the sum of autoconversion
and accretion processes. Speciﬁc formulations of each process are described in the following“Microphysical Processes” subsection.
The scheme utilizes a “dry” mixing ratio ( kg kg−1 ) to deﬁne the amount of water condensate.
For example, qc is the mass of cloud water per mass of dry air in the layer. The dry air density
ρ kg m−3 is calculated as ρ = P/(Rair T ), where P is the pressure in Pa, and the gas constant of
air Rair = 287.04 J kg−1 K−1 . A condensate mass is obtained by multiplying the mixing ratioby the air density. (e.g., the mass of ice mi = ρqi ). A number concentraion is in units m−3 .
Hereafter, unless stated otherwise, the cloud variables qc , qi , Nc , and Ni represent grid-averaged′
values; prime variables represent mean in-cloud quantities (e.g., such that qc = Cqc , where C
is cloud fraction). Note that qv ′ ̸= qv . qv ′ for ice clouds is determined as described in pmlsc
section. The sub-grid scale variability of water content within the cloudy area is not consideredat present.
3.5.2

Microphysical Processes

The time evolution of qi by microphysical processes is written in symbolic form as follows.

70



∂qi
∂t












∂qi
∂qi
∂qi
∂qi
+
+
+
∂t esnw∂t fallin∂t fallout∂t hom
micro








∂qi
∂qi
∂qi
∂qi
+
+
+
+
∂t het∂t dep∂t rim∂t mlt=(3.140), where t is time. The terms of the right hand side denote evaporation of snow (sunscriptesnw), ice fall in from above layers (subscript fallin), ice fall out to below layers (subscript
fallout), homogeneous nucleation (subscript hom), heterogeneous nucleation (subscript het),
deposition/sublimation (subscript dep), riming (subscript rim), and melting (subscript mlt).
Similarly, the time evolution of qc by microphysical processes is


∂qc
∂t




=
micro

∂qc
∂t




+

hom

∂qc
∂t




+

het

∂qc
∂t




+
rim



∂qc
∂t


+
evap

∂qc
∂t




+

auto

∂qc
∂t


(3.141)accr, where the terms on the right hand side are homogeneous nucleation, heterogeneous nucleation, riming, evaporation (subscript evap), autoconversion (subscript auto), and accretion(subscript accr). The formulations of these processes are detailed in the following subsections.
The conversion terms of all processes are calculated at every layer downward from the top
layer (k=kmax) to the bottom layer of the column (k=1). k is the vertical level increasing withheight, i.e., k+1 is the next vertical level above k.
The changes in the temperature of a layer is treated consistent with the phase-change ofwater.


∂T
∂t




=
phase change

∂T
∂t




+
vapor↔liquid

∂T
∂t




vapor↔solid

∂T
∂t



!

+


(3.142)liquid↔solid

with


∂T
∂t





∂T
∂t



∂T
∂t

Lv
=
cp
vapor↔liquid





Ls
=
cp
vapor↔solid

Lf=cpliquid↔solid



∂qi
∂t







∂qi
∂t


+

hom

∂qc
∂t


+
evap



∂qi
∂t

∂qr
∂t


+

esnw




+
het



∂qi
∂t

∂qi
∂t

(3.143)erain!
(3.144)

dep




+

rim

∂qi
∂t




,

(3.145)mlt

where Lv , Ls , and Lf is the latent heat of vaporization, sublimation, and fusion, respectively.
Cp is the speciﬁc heat of moist air at constant pressure.
3.5.2.1

Ice Properties

The formulation of the ice conversion terms requires parametrization of the mass, fall speed
and particle size distributions of ice. These are described ﬁrst and then subsequently used toderive the conversion terms.
71

The ice particle size distribution is parametrized as
Ni (D) = Ni0 exp(−0.1222(T − T0 )) exp (−Λi D) ,(3.146)

where D is the equivolume diameter of the particle in m, Ni0 = 2.0 × 106 m−4 , T is thetemperature in K, and T0 = 273.15 K. Λi represents the slope of the exponential distribution.
The temperature function exp(−0.1222(T − T0 )) represents the fact that ice particles tend tobe smaller at lower temperatures, and is an implicit way of parametrizing aggregation.
The mass of an ice particle is parametrized as a function of Dmi (D) = aDb(3.147)

where a = 0.069 kg m−2 and b = 2.0.
The fall-speed of an ice particle at an air density of ρ0 = 1 kg m−3 isvi (D, ρ0 ) = cDd(3.148)where c = 25.2m0.473 s−1 and d = 0.527.
At low air densities a particle will fall faster than at high air densities. Considering suchventilation eﬀect, the fall-speed of a particle at arbitrary air density ρ isvi (D, ρ) = (ρ0 /ρ)0.4 vi (D, ρ0 )(3.149)The combination of the size distribution, mass and velocity relationships yields a fall-speed
and ice water content relationship.
For a given ice content and temperature, Λi can be calculated by integrating (A.2) acrossthe particle size distribution (A.1). This gives the result that, for a given temperature, Λi is
proportional to the inverse cube root of the ice water content.
Λi =2aNi0 exp(−0.1222(T − T0 ))
mi

1
3

(3.150)3.5.2.2

Evaporation of Rain and Snow r
The evaporation rate of rain ∂q
is expressed as
∂t
erain



∂qr
∂t


=
erain

1
Fr
kE (q w − qv )
ρ∆z
VT r

(3.151), where Fr denotes the net accumulation of rain water at the layer in kg m−2 s−1 , VT r theterminal velocity, and kE the evaporation factor (VT r = 5 m s−1 and kE = 0.5). q w correcponds
to the saturation water vapor mixing ratio at the wet-bulb temperature. The evaporation occurs
only when q w − qv > 0.
 iSimilary to this, the evaporation rate of falling ice ∂q
is expressed as
∂t
esnw



∂qi
∂t


= kE (q w − qv )
esnw

Fi
VT r

(3.152)where Fi denotes sedimentation of cloud ice from above layers. VT s is set to 5 m s−1 .
72

3.5.2.3

Ice Fall

The total ice ﬂux from the layer ‘k’ isZ ∞Fi |k =Ni (D)mi (D)vi (D)dD.

(3.153)0

The fraction of ice ﬂux from level the ‘k’ to the below level ‘kk’ (1 <= kk < k) iceweight|k,kk ,is given as
R f (zm(k)−zm(kk))
0

R f (zm(k)−zm(kk+1))
Ni (D)mi (D)vi (D)dD − 0
Ni (D)mi (D)vi (D)dD
R∞,
0 Ni (D)mi (D)vi (D)dD(3.154)

where zm(k) is the middle of the height of the layer k, and f(dz) is the ice size which falls
the distance dz in a single time step.
The net ice fall out from the layer is


∂qi
∂t


=−
fallout∆tFiρ∆z(3.155).
The net ice fall in to the layer ‘k’ is

3.5.2.4∂qi∂tl=kmax∆t X=
Fi |k=l × iceweight|l,kρ∆zfallin(3.156)l=k+1

Homogeneous Nucleation

This term simply converts all liquid cloud water to ice if the temperature is less than a given
threshold of 233.15 K.


3.5.2.5

∂qi
∂t




=−

hom

∂qc
∂t


=
hom

qc
∆t

(3.157)Heterogeneous Nucleation

A Spectral Radiation-Transport Model for Aerosol Species (SPRINTARS; Takemura etal. 2000, 2002, 2005, 2009) coupled with MIROC6 explicitly predicts the masses and number
concentrations for aerosol species. Heterogeneous freezing of cloud droplets takes place throughcontact and immersion freezing on ice cucleating particles (INPs), which are parameterized according to Lohmann and Diehl (2006) and Diehl et al. (2006). Soil dust and black carbon are
assumed to act as INPs. Ratios of activated INPs to the total number concentration of soil
dust and black carbon for the contact freezing and the immersion/condensation freezing are
based on Fig. 1 in Lohmann and Diehl (2006). Using the number of INPs (Nnuc ) predicted in
SPRINTARS, the rate of heterogeneous freezing is diagnosed as follows.


∂qi
∂t




=−

het

∂qc
∂t


= max{Nnuc Wnuc0 ,het

The weight of nucleated drop, Wnuc0 , is set to 1.0 × 10−12 .
73

qc
}
∆t

(3.158)

3.5.2.6

Deposition/Sublimation

A single ice particle grows or disappears by water vapor diﬀusion according to the following
equation:∂mi (D)= {4πC (Si − 1) F } / [{Ls /(Rv T ) − 1} Ls / (ka T ) + Rv T / (XPsati )]∂t(3.159)i (D)where ∂m∂t
is the rate of change of the particle mass; (Si − 1) is the supersaturation of
the atmosphere with respect to ice; Rv is the gas constant for water vapour; ka is the thermal
conductivity of air at temperature T, X is the diﬀusivity of water vapour; Psati is the saturated
vapour pressure over ice; Ls is the latent heat of sublimation of ice; C is a capacitance term andF is a ventilation coeﬃcient. C is assumed to appropriate to spheres, so is equal to D/2 . F isgiven by Pruppacher and Klett (1997) as F = 0.65 + 0.44Sc1/3 Re1/2 , where Sc is the Schmidt
number, equal to 0.6, and Re is the Reynolds number, v(D)ρD/µ, where µ is the dynamic
viscosity of air.
 i
is obtained as
Integrating ice size distribution, ∂q
∂t

dep



∂qi
∂t



1
=
ρ
dep

Z

∂mi (D)
N (D)dD∂t(3.160)The ice grows or disappears depending on the sign of (Si − 1).
1. (Si − 1) > 0
The ice grows (deposition). If qc exists in the same grid, qc is evaporated as fast as thedeposition process (Wegener–Bergeron–Findeisen process).


∂qc
∂t




=−

evap

∂qi
∂t


(3.161)dep

The basis of this theory is the fact that the saturation vapor pressure of water vapor withrespect to ice is less than that with respect to liquid water at the same temperature. Thus,
within a mixture of these particles, the ice would gain mass by vapor deposition at the expenseof the liquid drops that would lose mass by evaporation.
2. (Si − 1) < 0
The ice disappears (sublimation).
3.5.2.7

Cloud water Collection by Ice (Riming)

Riming process (the ice crystals settling through a population of supercooled cloud droplets,
freezing them upon collision) is based on the geometric sweep-out integrated over all ice sizes(Lohmann 2004):


∂qi
∂t




=−

rim

∂qc
∂t


=

πESW n0S aqc Γ(3 + b)
(3+b)

4λS

rim

74



ρ0
ρ

0.5
(3.162)

where n0S = 3 × 106 m−4 is the intercept parameter, λS is the slope of the exponential
Marshall-Palmer ice crystal size distribution, a = 4.84, b = 0.25, and ρ0 = 1.3 kg m−3 is thereference density. Γ is the gamma fanction. The collection eﬃciency Esw is highly dependent
on the cloud droplet and snow crystal size (Pruppacher and Klett 1997). The size-dependent
collection eﬃciency for aggregates is introduced as obtained from laboratory results by Lew etal. (1986) (simulation ESWagg).
agg
ESW= 0.939St2.657(3.163)The Stokes number (St) is given bySt =2 (Vt − vt ) vt
Dg

(3.164)

Vt is the snow crystal terminal velocity, and D is the maximum dimension of the snow crystal.
vt is the cloud droplet terminal velocity. g is the acceleration due to gravity.
3.5.2.8

Ice Melt

Since this term is essentially a diﬀusion term, although of heat instead of moisture, its form
is very similar to that of the deposition and evaporation of ice term. The rate of change of ice
mass of a melting particle is given by:


∂qr
∂t




=−
mlt

∂qi
∂t


= 4πCF {ka /Lm (T w − T0 )}(3.165)mlt,where Lm is the latent heat of melting of ice, T w is the wet-bulb temperature of the air andT0 = 273.15K is the freezing point of ice. Ice melt occurs when T w − T0 > 0. The capacitanceterm, C, is considered to be that for spherical particles. Hence C = D/2. The ventilation factor,
F is considered to be the same as in the deposition/sublimation process.
3.5.2.9

Warm Rain Cloud Microphysics

We assume Nc is the number of aerosols activated as droplets. The nucleation of cloud
droplets is predicted in the aerosol module SPRINTARS (Takemura et al. 2000; 2002; 2005;2009) based on the parameterization by Abdul-Razzak and Ghan (2000), which depends on
the aerosol particle number concentrations, size distributions and chemical properties of of eachaerosol species, and the updraft velocity.
The autoconversion term following Berry (1968) is a function of qc and Nc .


∂qr
∂t




=−

auto

∂qc
∂t


=
auto

1 b1 × m′2cNcρ b2 + b3 m′

(3.166)c

The parameters are set as b1 = 0.035, b2 = 0.12, b3 = 1.0 × 10−12 . The eﬀect of aerosol-cloudinteraction on cloud lifetime is taken into account by the dependency on Nc .
The accretion term is given as


∂qr
∂t




=−

auto

∂qc
∂t


=
auto

1
qc qr
ρ

(3.167)

Rain water qr into the layer is from above the layer. qr is treated as a diagnostic variables:
qr falls out to surface within the time step.
75

3.5.2.10

Total Precipitation

The total amount of precipitation at a certain pressure level, p, is obtained by integrating
the relevant processes from the top of the model (p = 0) to the respective pressure level. The
ﬂuxes of rain and ice kgm−2 s−1 are then expressed as

1
Prain (p) =g

1
Pice (p) =g

Z p 
0



Z p
0

∂qi
∂t

∂qr
∂t





+
auto


−

fallin



∂qi
∂t

∂qr
∂t




+

accr




+
fallout

76

∂qi
∂t

∂qr
∂t





−
mlt


−

rim



∂qr
∂t

∂qr
∂t



dp

(3.168)

revap


−

mlt

!



∂qr
∂t




dp(3.169)

esnw

3.6
3.6.1

Radiation Scheme
Summary of the Radiation Flux Calculation

The radiation scheme in the MIROC was created based on the Discrete Ordinate Method and thek-distribution Method (Nakajima et al., 2000), and updated by Sekiguchi and Nakajima (2008).
The scheme calculates the value of the radiation ﬂux at each level by considering the absorption,
emission, and scattering processes of terrestrial and solar radiation by gases and clouds/aerosols.
The main input data are temperature T , speciﬁc humidity q, cloud water l, and cloud cover C.
The output data are shortwave or longwave upward and downward radiation ﬂuxes F ∓ , and+derivative coeﬃcient to surface temperature dF ∓ /dTg , surface downward radiation ﬂux Fsf,
vis
and 0.5 and 0.67 µm optical thickness τ .
The calculation is separated for several wavelength bands. It is further divided into several
sub-channels, based on the k-distribution method. As for gaseous absorption, the line absorptionin H2 O, CO2 , O2 , O3 , N2 O, CH4 , the continuous absorption in H2 O, CO2 , O2 , O3 , and the CFC
absorption are incorporated. As for scattering, Rayleigh scattering of gases and scattering by
cloud and aerosol particles are considered.
Major subroutines used to calculate the radiation ﬂux in SUBROUTINE:[DTRN31] of pradt.F
are as follows.
1. Calculate the Planck function from atmospheric temperature SUBROUTINE:[PLANKS, PLANKF]
2. Calculate the optical thickness to the gas in each sub-channel SUBROUTINE:[PTFIT2]
3. Calculate the optical thickness to the CFC absorption SUBROUTINE:[CNTCFC2]
4. Calculate the optical thickness to aerosol, Rayleigh scattering, and cloud SUBROUTINE:[SCATAE,SCATRY, SCATCL]
5. Expand the Planck function by optical thickness for each sub-channel SUBROUTINE:[PLKEXP]
6. Calculate the transmission coeﬃcient (T), reﬂection coeﬃcient (R) and source function(S) SUBROUTINE:[TWST]
7. Make T, R, and S matrixes for maximal/random approximation SUBROUTINE:[RTSMR]
8. Calculate the radiation ﬂux by adding method SUBROUTINE:[ADDMR, ADDING]
To account for the partial coverage of clouds, the transmission and reﬂection coeﬃcients andsource functions for each layer are calculated at weighted average of the cloud cover, separately
for cloud cover and clear-sky conditions. The cloud cover of the cumulus is also considered. In
addition, it also performs several adding and calculates the clear-sky radiation ﬂux.
3.6.2

Wavelength and Sub-Channel

The basics of radiative ﬂux calculations are represented by Beer-Lambert’s Law.
F λ (z) = F λ (0)exp(−k λ z)

(3.170)F λ is the radiant ﬂux density at the wavelength of λ and k λ is the absorption coeﬃcient. In
order to calculate the radiative ﬂuxes related to the heating rate, the integration operation withrespect to the wavelength is required.
Z
F (z) =ZλF (z)dλ =F λ (0)exp(−k λ z)dλ(3.171)

However, it is not easy to calculate this integration precisely because the absorption andemission of radiation by gas molecules have the complicated wavelength dependence of the77

Figure 4: Flowchart of SUBROUTINE:[DTRN31]
absorption line attributed to the structure of the molecule. The k-distribution method is amethod designed to make the relatively precise calculation easier. Within a certain wavelength
range, considering the density function F (k) for λ of the absorption coeﬃcient of k, the above
formula is approximated as follows,Z

Z
F (0)exp(−k z)dλ ≃λλF̄ k (0)exp(−kz)F (k)dk(3.172)where F̄ k (0) is the ﬂux averaged over a wavelength having the absorption coeﬃcient in thiswavelength k in z = 0.
If F̄ k (0) and F (k) are a relatively smooth functions to the k,ZF λ (0)exp(−k λ z)dλ ≃XF̄ i (0)exp(−k i z)F i(3.173)the formula, as such above, can be relatively precisely calculated by the addition of a ﬁnite
number (sub-channels) of exponential terms. This method has furthermore the advantage easy
to consider the absorption and scattering at the same time.
In the MIROC 6.0, by changing the radiation parameter data, the calculations can be performed at various wavelengths. In the standard version, the wavelength range is divided into 29
parts. In addition, each wavelength range is divided into 1 to 6 sub-channels (corresponding to
the i in the above formula). There are 111 channels in total. The wavelength range is divided
by the wavenumber ( cm−1 ), 1, 250, 400, 530, 610, 670, 750, 820, 980, 1175, 1225, 1325, 1400,
2000, 2500, 3300, 3800, 4700, 5200, 6000, 10000, 12750, 13250, 14750, 23000, 30000, 33500,
36000, 43500, 50000. Additionally, a chemical version is also with 37 bands and 126 channels
78

for chemical transport model and the boundary of the shortwave region is also changed to 54000
cm−1 .
3.6.3

Calculation of the Planck Function

In this section, SUBROUTINE:[PLANKS, PLANKF] in pradt.F is described.
The Planck function B̄ w (T ), integrated in each wavelength range, is evaluated by the following formula.

B̄ w (T ) = λ−2 T exp( 5
X

Bnw λ̄w T

−n)(3.174)n=1

where λ̄w is the averaged wavelength of the wavelength range, Bnw is the parameter determined by function ﬁtting. This is calculated to the atmospheric temperature of each layer Tl ,and the boundary atmospheric temperature of each layer Tl+1/2 , surface temperature Tg andtemperature 1K higher than surface temperature Tg+1K . The calculations are performed foreach wavelength and each layer. In the following description, the subscript of the wavelength
range w is omitted.
3.6.4

Calculation of the Optical Thickness to Gas Absorption

In this section, SUBROUTINE:[PTFIT2] in pradt.F is described.
The optical thickness of the gas absorption (the line and continuum absorption are uniﬁed)
τ KD is expressed as follows by using the index m as the type of molecules.
τ KD =X

k (m) C (m)(3.175)m=1where k (m) is the absorption coeﬃcient of the molecule m, which is diﬀerent for each subchannel and determined as a function of temperature T and atmospheric pressure p. C (m)represents the amount of gas in the layer represented by mol/cm2 /km, calculated by using thegas concentration r(m) in ppmv ( C (m) = 10−1 r(m) ρdz ). In the MIROC 6.0, the number ofthe considered molecule types m is 6 (1:H2 O, 2:CO2 , 3:O3 , 4:N2 O, 5:CH4 , 6:O2 ). Also, k (m) isrepresented as follows (the details are in Sekiguchi and Nakajima, 2008).
(m)
k (m) = exp log 10k2 + (A + BT ) log (T /Tref2 )




(m)
(m)
(m)
(m)
log 10 k1 − k2
log 10 k3 − k2
 / (Tref3 − Tref1 )




−
B=
T 3
Tref1
log Tref
log
Tref2
ref 2

(3.176)



A=



(m)
(m)
log 10 k3 − k2log (Tref3 /Tref2 )− BTref3(3.177)(3.178)(m)

Tref 1−3 are the reference temperatures prepared in advance (200, 260, 320 K), and k1−3 are
the absorption coeﬃcients when the reference temperatures Tref 1−3 is used (also ﬁtted at 26
atmospheric pressure grids).
79

When considering the absorption of H2 O, we calculate the optical thickness of the selfbroadening and add τ self .
τ KD(H2 O) = τ KD(H2 O) + τ self

(3.179)2

τk (H2 O self) C (H2 O)= (H O)
C 2 + ρdz105

self(3.180)k (H2 O self) is calculated in the same way as k (m) . The self-broadening absorption coeﬃcientsin the reference temperatures Tref 1−3 are prescribed and dependent on the pressure. In theabove formula, 105 is multiplied to convert the unit from km to cm. This calculation is done foreach sub-channel and each layer.
3.6.5

Calculation of the Optical Thickness to CFC Absorption

In this section, SUBROUTINE:[CNTCFC2] in pradt.F is described.
The optical thickness of the CFC absorption τ CF C is considered for several types of CFCsm.
τ CF C =

X

10k

(m)

r(m) ρ∆z10−1(3.181)m

In MIROC 6.0, the number of the considered CFCs m is 28 (1:CFC-11, 2:CFC-12, 3:CFC-13,
4:CFC-14, 5:CFC-113, 6:CFC-114, 7:CFC-115, 8:HCFC-21, 9:HCFC-22, 10:HCFC-123, 11:HCFC-124,
12:HCFC-141b, 13:HCFC-142b, 14:HCFC-225ca, 15:HCFC-225cb, 16:HFC-32, 17:HFC-125, 18:HFC-134,
19:HFC-134a, 20:HFC-143a, 21:HFC-152a, 22:SF6 , 23:ClONO2 , 24:CCl4 , 25:N2 O5 , 26:C2 F6 ,27:HNO4 , 28:SF5 CF3 ). In the above formula, 10−1 is multiplied to convert from km to cm,
and from ppmv to ratio. This calculation is done for each sub-channel and each layer. This
calculation is performed for each layer and the wavelength range from about 540 to 1800 cm−1 .
3.6.6

Optical Thickness to Scattering and Scattering Moment

Calculate the optical thickness of scattering and the scattering moment. These calculations are
performed for each wavelength and each layer. The optical parameters for the particle matter(p)
qm are prepared, including the extinction coeﬃcient (m = 1) including the scattering andabsorption process and the absorption coeﬃcient (m = 2) the moments of the volume scattering
phase function (m = 3-4: ﬁrst-second order).
3.6.6.1

Aerosol

In this section, SUBROUTINE:[SCATAE] in pradt.F is described.
ae , the scatThe optical thickness τ ae , the part of the optical thickness due to absorption τabaetering moment Qm for aerosol areτ ae =X (p)q1,n r(p) × ρ∆z10−1p80(3.182)

ae
τab
=

X (p)
q2,n r(p) × ρ∆z10−1(3.183)

p

Qae
m =

X

(p) (p)qm,nr × ρ∆z10−1 ( m ≥ 3)

(3.184)p

p is the aerosol type, and r(p) is volume mixing ratio of the particle. The optical parameters(p)for the particle qm,n depend on the mode radius index n prescribed for each particle (IRA).
In the MIROC 6.0, the number of the considered aerosol types p 15 (1-6:soil dust (bin1-6),
7:carbonaceous (BC/OC=0.3), 8:carbonaceous (BC/OC=0.15), 9:carbonaceous (BC/OC=0),
10:black carbon (external mixture), 11:sulfate, 12-15:sea salt (bin 1-4)).
If the aerosol radius is used, the optical thickness τ ae , the part of the optical thickness due toae , and the scattering moment Qae for the hygroscopic aerosols (e.g., carbonaceous,absorption τabmsulfate, sea salt) areτ ae =

i
Xh
(p)
(p)
(1 − F Xae ) q1,nf it r(p) + F Xae q1,nf it+1 r(p) × ρ∆z10−1(3.185)

p

ae
=
τab

i
Xh
(p)
(p)
(1 − F Xae ) q2,nf it r(p) + F Xae q2,nf it+1 r(p) × ρ∆z10−1(3.186)

p

Qae
m =

Xh

i
(p)
(p)
(1 − F Xae ) qm,nf it r(p) + F Xae qm,nf it+1 r(p) × ρ∆z10−1 ( m ≥ 2)

(3.187)p(ref )F Xae = RH − RHnf it







1
(ref )(ref )
RHnf it+1 − RHnf it(3.188)(ref )

where RH is the local relative humidity and RHnf it is the relative humidity given in theparameter and nf it is the number of the prescribed relative humidity closest to the RH. nf it
and F Xae are calculated in the SUBROUTINE:[RMDIDX] in pradt.F and determined in advance.
In the above formulas, 10−1 is multiplied to convert from km to cm, and from ppmv to ratio.
3.6.6.2

Rayleigh Scattering

In this section, SUBROUTINE:[SCATRY] in pradt.F is described.
The optical thickness τ r of Rayleigh scattering and the part of the optical thickness due tor areabsorption τabτr =er qmol1 dp
pST D(3.189)r
τab=er qmol2 dp
pST D(3.190)81

pST D = 1013.25(3.191)

where er is the Rayleigh scattering coeﬃcient, qmolm is the moments of the phase function.
These calculations are performed up to m = 2. Also, this is added to the optical thickness forthe aerosol.

3.6.6.3

τ ae+r = τ ae + τ r(3.192)

ae+r
ae
r
τab
= τab
+ τab

(3.193)Cloud

In this section, SUBROUTINE:[SCATCL] in pradt.F is described.
cl , and the
The optical thickness τ cl , the part of the optical thickness due to absorption τabclscattering moment Qm for cloud areτ cl =X (ct)
q1,n r(ct) × ρ∆z10−1(3.194)

ct

cl
=
τab

X (ct)
q2,n r(ct) × ρ∆z10−1(3.195)

ct

Qcl
m =

X

(ct)
qm,n
× r(ct) ρ∆z10−1 ( m ≥ 3)

(3.196)ct

ct is the cloud particle type (1:liquid cloud, 2:ice cloud). The optical parameters for the(ct)particle qm,n depend on the mode radius index n prescribed for each particle (IRC). If the cloud
cl ,radius is used, the optical thickness τ cl , the part of the optical thickness due to absorption τaband the scattering moment Qcl
m for cloud areτ cl =

i
Xh
(ct)
(ct)
(1 − F Xcl ) q1,nf it r(ct) + F Xcl q1,nf it+1 r(ct) × ρ∆z10−1(3.197)

ct

cl
τab
=

i
Xh
(ct)
(ct)
(1 − F Xcl ) q2,nf it r(ct) + F Xcl q2,nf it+1 r(ct) × ρ∆z10−1(3.198)

ct

Qcl
m =

Xh

i
(ct)
(ct)
(1 − F Xcl ) qm,nf it r(ct) + F Xcl qm,nf it+1 r(ct) × ρ∆z10−1 ( m ≥ 3)

(3.199)ct

(ref )F Xcl = R(ct) − Rnf it






82


1
(ref )(ref )
Rnf it+1 − Rnf it(3.200)(ref )

where R(ct) is the calculated mode radius and Rnf it is the mode radius given in the parameterand nf it is the number of the prescribed mode radius closest to the R(ct) . nf it and F Xcl arecalculated in the subroutine SUBROUTINE:[RMDIDX] in pradt.F and determined in advance. In
the above formulas, 10−1 is multiplied to convert from km to cm, and from ppmv to ratio.
Finally, the total optical thickness for particle scattering, Rayleigh scattering and absorption
τ p and the contribution of scattering τ scat are obtained as follows.
τ P = τ cl + τ ae+r

(3.201)



ae+r
cl
τ scat = τ P − Tab+ Tab(3.202)

In addition, the moments of the normalized phase function G are calculated up to the threeorders. The zeroth moment G1 is trivial from the normalization condition of the phase function.
The ﬁrst and second moments G2 , G3 , are referred as the asymmetry factor g and the truncation
factor f .
G1 = 1.0

Gm−1 =aeQcl
m + Qm(m ≥ 3),
τ scat(3.203)G2 = g,

G3 = f(3.204)This calculation is divided into the cloudy, clear sky and cumulus conditions. In the cloudy
and cumulus conditions, τ cl in the 0.5 and 0.67 µm regions is as recorded as τ vis in subroutine
DTRN31.
**Rct is calculated in SUBROUTINE:[RADFLX] as follows.

R(ct) =

3 ρr(ct)(ct)
4π ρ(ct)
w nc!1/3(3.205)(ct)

ρw is the liquid or ice density. rct is the amount of the liquid or ice cloud and calculated
as follows.
(ct)r(ct) =

(ct)

Cst rst + Ccu rcu
1 − (1 − Cst ) (1 − Ccu )(3.206)

C is the area of the cloud, and the subscript st and cu mean the stratus and cumulus. When(ct)
nc is the number density ofcloud particles.
(ct)
rst,cu is the small amount in the stratosphere, it is reset to 0.

n(liq)= maxc
n(ice)
= max
c

liq
qae
pNA
(liq)
, fliq nmin
RTv (18 × 10−3 Rv /R)!

ice pNqae(ice)A, (1 − fliq ) nmin
RTv (18 × 10−3 Rv Rv /R)

83(3.207)(3.208)liq
where qaeis the mixing ratio of the aerosol particles calculated by the SPRINTERS and(ct)
converted to the number concentration, and nmin is the minimum number of the cloud particles.
(ct)
and fliq is liquid fraction. Also, nc is calculated as follows when using OPT AECL SIMPLE.

(ct)

n(ct)
=
c

εna nmax(ct)(3.209)εna + nmax

where na is the number density of aerosol particles give as an external condition, and ε and(ct)
nmax are constants. fliq is calculated by the following formula using the amount of cloud waterw (0 ≤ fliq ≤ 1).

fliq =
3.6.7

wst fliq,st + wcu fliq,cu
wst + wcu

(3.210)Total Optical Thickness

All optical thickness including gaseous band absorption, and scattering is,τ = τ KD + τ CON + τ P(3.211)

where because τ KD is diﬀerent for each subchannel, the calculation is done for each subchannel and each layer, and divided into the cloudy, clear sky, and cumulus conditions.
3.6.8

Expansion of the Plank Function

In this section, SUBROUTINE:[PLKEXP] in pradt.F is described.
In each layer, the Planck function B is expanded as follows and the expansion coeﬃcientsb0 , b1 , and b2 , are obtained.
2
B τ ′ = b0 + b1 τ ′ + b2 τ ′(3.212)

Here, as B (τ ′ ), B at the top of each layer (the boundary with the top layer) is used, and asB(τ ), B at the bottom edge of each layer (the boundary with the layer below), and as B(τ /2),
the B at the representative level of each layer.

b0 = B(0)(3.213)b1 = (4 B(τ /2) − B(τ ) − 3 B(0))/τ

(3.214)b2 = 2( B(τ ) − B(0) − 2 B(τ /2))/τ 2(3.215)This calculation is done for each sub-channel and each layer and divided into the cloudy,
clear sky and cumulus conditions.

84

Figure 5: Second-order expansion using the optical thickness of the plank function
3.6.9

Transmission and Reﬂection Coeﬃcients, and Source Function

In this section, SUBROUTINE:[TWST] in pradt.F is described.
Using the obtained optical thickness τ , optical thickness of scattering τ scat , scattering moments g, f , expansion coeﬃcients for Planck function bn , and solar incidence angle factor µ0 , thetransmission coeﬃcient T , reﬂection coeﬃcient R, downward radiation source function ϵ+ , andthe upward radiation source function ϵ− are founded, by assuming a uniform layer and using
the two-stream approximation.
The single-scattering albedo ω is,ω = τ scat /τ

(3.216)The optical thickness τ ∗ , the single-scattering albedo ω ∗ , and asymmetric factor g ∗ , correctedby the contribution from the forward scattering factor f are,τ ∗ = τ (1 − ωf )(3.217)(1 − f )ω
1 − ωf(3.218)g−f
1−f(3.219)ω∗ =g∗ =µ is a two-stream directional cosine.

µ≡

1
1
√ ,
3 1.66(shortwave, longwave)W − ≡ µ−1/285(3.220)(3.221)

Furthermore,P̂ ± = ω ∗ W −2 1 ± 3g ∗ µ2 /2(3.222)Ŝs± = ω ∗ W − (1 ± 3g ∗ µ0 µ)(3.223)can be determined as above as a normalized scattering phase function.

X=µ

−1



− P̂

+

− P̂

+


(3.224)= µ−1 − 3ω ∗ W −2 µ2 g ∗
Y = µ−1 − P̂ + + P̂ +(3.225)= µ−1 − ω ∗ W −2σ̂S+ = ŜS+ + ŜS−

(3.226)= ω∗W −σ̂S− = Ŝs+ − ŜS−

(3.227)= 3ω ∗ µ0 W − µg ∗Using the above formula, the reﬂectance R and transmission T become
(
) (
)
∗
∗
X 1+e−λτ −λ 1−e−λτ(3.228)AA = X 1+e−λτ ∗ +λ 1−e−λτ ∗
(
) (
)

(
) (
)
∗
∗
X 1−e−λτ −λ 1+e−λτ(3.229)BB = X 1−e−λτ ∗ +λ 1+e−λτ ∗
(
) (
)

λ=

√

(3.230)XY

R = 12 (AA + BB)(3.231)

T = 12 (AA − BB)(3.232)

Next, we ﬁnd the source function derived from the Planck function.

∗b̂n = 2π (1 − ω ) W

−

86



1
1 − ωfnbn(3.233)

The expansion coeﬃcients of the radiation source function can be found from the above
formulas.
D2± = b̂Y2

(3.234)2b̂2
D1± = b̂Y1 ∓ XY

(3.235)2b̂2b̂1
D0± = b̂Y0 + XY
2 ∓ XY(3.236)

D± (τ ∗ ) = D0+ + D1+ τ ∗ + D2+ τ ∗2

(3.237)

The radiation source function derived from the Planck function ϵ̂±
A is
−
−
−
∗
ϵ̂−
A = D0 − RD0 − T D (τ )(3.238)

+
−
∗
+
∗
ϵ̂+
A = D (τ ) − T D0 − RD (τ )(3.239)On the other hand, the radiation source function of the solar-induced radiation isϵ̂+
S = Fsol



τ∗
τ∗
−µ
−µ
+
+
−
0
0
Vs e
− T Vs − RVs e(3.240)ϵ̂+S = Fsol



τ∗
τ∗
−µ
−µ
+
−
+
0
0
Vs e
− T Vs − RVs e(3.241)

Here, Qγ and Vs± are expressed by the following formulas, and Fsol is solar irradiance.
h

− i
σ̂S
+
Vs± = 12 Qγ ± Qγ
µ0
X

Qγ =

−
+
X σ̂S
µ0 +µ−10 σ̂S
XY µ0 +µ−10(3.242)(3.243)The direct solar transmission is also calculated in this subroutine.
∗

Exdir = e−τ /µ0(3.244)This calculation is done for each sub-channel and each layer and divided into the cloudy,
clear sky, and cumulus conditions.

873.6.10

T, R, S Matrixes for Maximal/Random Approximation

In this section, SUBROUTINE:[RTSMR] in pradt.F is described.
In this subroutine, T, R, S matrixes for maximal/random approximation is made. The
radiation source function, which is the sum of both the plank function and the solar incident
origins, is−(cloud)ϵ−(cloud) = ϵ̂S−(clear)ϵ−(clear) = ϵ̂S−(clear)T r(clear) + ϵ̂A+(cloud)ϵ+(cloud) = ϵ̂S+(clear)ϵ+(clear) = ϵ̂S−(cloud)

T r(cloud) + ϵ̂AC(3.245)(1 − C)(3.246)−(cloud)

T r(cloud) + ϵ̂A−(clear)T r(clear) + ϵ̂AC(3.247)(1 − C)(3.248)T r is the direct solar transmission for maximal/random approximation and calculated asfollows in SUBROUTINE:[DTRN31].
(cloud) (3)(clear)(1)
Tr(cloud)
=
Ex
B
+
Ex
1
−
B
n
n
n
n
n

(cloud)= T r(cloud)Exdir(cloud)
n
n

Exn+1



T rn(clear) = Ex(cloud)
1 − Bn(3) + Ex(clear)
Bn(1)
n
n

(clear)Exn+1(3.249)(3.250)(3.251)= T r(clear)
Exdir(clear)nn(3.252)(1−4)

is the cloud cover interaction index

Ex is the cumulative direct solar transmission. Bn
and calculated in SUBROUTINE:[BCVR]in pradt.F.
(1)(3.253)(2)(3.254)n−1 ,Cn )Bn = 1−max(C1−Cn−1n+1 ,Cn )Bn = 1−max(C1−Cn+1(3)(3.255)(4)(3.256)n−1 ,Cn )Bn = min(CCn−1

n+1 ,Cn )Bn = min(C
Cn+1

88

Next, T matrixes for maximal/random approximation are calculated.
T +(cloud,1) = T (cloud) B (3)(3.257)

T +(cloud,2) = T (cloud) 1 − B (1)T +(clear,1) = T (clear) 1 − B (3)(3.258)(3.259)T +(clear,2) = T (clear) B (1)(3.260)T −(cloud,1) = T (cloud) B (4)(3.261)T −(cloud,2) = T (cloud) 1 − B (2)T −(clear,1) = T (clear) 1 − B (4)

T −(clear,2) = T (clear) B (2)(3.262)(3.263)(3.264)

Also, R matrixes for maximal/random approximation are calculated.
R+(cloud,1) = R(cloud) B (3)(3.265)R+(cloud,2) = R(cloud) 1 − B (1)R+(clear,1) = R(clear) 1 − B (3)(3.266)(3.267)R+(clear,2) = R(clear) B (1)(3.268)R−(cloud,1) = R(cloud) B (4)(3.269)R−(cloud,2) = R(cloud) 1 − B (2)R−(clear,1) = R(clear) 1 − B (4)R−(clear,2) = R(clear) B (2)This calculation is done for each sub-channel and each layer.
89(3.270)(3.271)(3.272)

Figure 6: Schematic illustration of the adding method
3.6.11

Adding of Source Functions for Each Layer

In this section, SUBROUTINE:[ADDMR] and SUBROTINE:[ADDING] in pradt.F is described.
By using transmission coeﬃcient T , reﬂection coeﬃcient R, and radiation source function εin all layers, the radiation ﬂuxes u at each layer boundary can be obtained by using the adding
method. This means that the when two layers of T , R, ε are known, the T , R, ε of the wholecombined layer of the two layers can be easily calculated.
3.6.11.1

SUBROUTINE:[ADDMR]

In this subroutine, the maximal/random ﬂux in cloudy conditions is calculated by the adding
method and the T, R, and S matrixes are used for calculations.
First, the upward radiation source function the bottom layer is calculated.
In the shortwave region,

ϵN−(cloud)= W + αs µ0−(clear)= W + α s µ0

ϵN

 
1
µ

 
1
µ

90

−⟨τ ∗ ⟩/µ0 (cloud)F0 eN−⟨τ ∗ ⟩/µ0 (clear)

F0 eN(3.273)(3.274)⟨τ ∗ ⟩ indicates the total optical thickness τ ∗ of from the upper end of the atmosphere to the up∗
per end of the layer currently being considered and e−⟨τ ⟩/µ0 is calculated in SUBROUTINE:[ADDMR]
of pradt.F.
In the longwave region,−(cloud)ϵN−(clear)ϵN−(cloud)= ϵN+ W + 2π (1 − αs ) BN CN(3.275)−(clear)+ W + 2π (1 − αs ) BN (1 − CN )(3.276)

W + ≡ µ1/2

(3.277)

= ϵN

Here,

−
The reﬂectance R1,nand source function ϵ+1,n regarded from the ﬁrst to the n layers as asingle layer are


−1 

+
−
+
+
+ −
ϵ−
R1,n−1ϵ+n + ϵ1,n−1
1,n = ϵn + Tn 1 − Rn R1,n−1(3.278)


−1
−
−
−
Tn−
= Rn− + Tn+ 1 − Rn+ R1,n−1
R1,n−1R1,n(3.279)The upward and downward ﬂuxes at the bottom of the atmosphere are

−1 

−
−
+
+
−
ϵ
+
R
ϵ
R
=
1
−
R
u+
1,N −1 N
1,N −1
1,N −1 N
N +1/2(3.280)


−1 

+ +
−
+ −
ϵ
+
R
ϵ
R
u−
=
1
−
R
N
N 1,N −1
N 1,N −1
N +1/2(3.281)

Also, upward and downward ﬂuxes at the boundary between layers n and n + 1 are


−1 
−
−
−
− u−
− + ϵ++(3.282)
R
R
T
+
R
ϵ
u+
=
1
−
R
1,n
1,n−1 n1,n−1 n n+1/2
1,n−1 n
n+1/2

−1 

−
+ R−
− u−
+ ϵ+
u−
=
1
−
R
T
+
R
+
ϵ
n 1,n−1
n n+1/2
n 1,n−1nn+1/2(3.283)

However, the upward and downward ﬂux at the upper end of the atmosphere is as follows.
u+1/2 = 0(3.284)− −−u−
1/2 = T1 u1+1/2 + ϵ1

(3.285)

Finally, since this ﬂux is scaled, we rescaled and added the direct solar incidence to the ﬁnd
the radiation ﬂux. Furthermore, the ﬂux in the cloudy area and the clear sky area is added.
91

+

+
Fn+1/2= WW+(cloud)

un+1/2+(clear)+ un+1/2−
Fn+1/2= WW+−(cloud)

un+1/2−⟨τ ∗ ⟩/µ (cloud)+ µ0 F0 en+1/2 0−(clear)(3.286)

+ un+1/2(3.287)+
Also, surface direct downward radiation ﬂux Fsfis−⟨τ ∗ ⟩/µ0 (cloud)−⟨τ ∗ ⟩/µ0 (clear)+
Fsf= µ0 F0 eN+ eN(3.288)This calculation is done for each sub-channel.
3.6.11.2

SUBROUTINE:[ADDING]

Since the maximal/random approximation cannot be used under the clear sky condition,
this subroutine is used to calculate the ﬂux.
First, the radiation source function, which is the sum of both the plank function origin andthe solar incident origin, is∗−⟨τ ⟩/µ0+ ϵ̂±ϵ± = ϵ̂±ASe(3.289)There are layers 1, 2, . . . , N from the top. The surface layer is considered to be a single layer
and the N layer. Given the reﬂectance and source function of the layers from the n to the Nlayer as a single layer Rn,N , ϵ−n,N ,Rn,N = Rn,N + Tn (1 − Rn+1,N Rn )−1 Rn+1 Tn(3.290)−1−−
Rn+1,N ϵ+ϵ−
n + ϵn,N
n,N = ϵn + Tn (1 − Rn+1,N Rn )(3.291)

Rn,N = Rn,N + Tn (1 − Rn+1,N Rn )−1 Rn+1 Tn(3.292)−1− + T (1 − R+ + ϵ−
ϵ−
=
ϵ
R
)
R
ϵ
n
n
n+1,N
n+1,N
n
n
n,N
n,N

(3.293)This can be solved by n = N − 1, . . . , 1 in sequence, starting from the values at the surface
RN,N , ϵ−N,N .
RN,N = RN = αs(3.294) ∗ϵN,N = W + αs µ0 µ1 e−⟨τ ⟩/µ0 F0 + 2π (1 − αs ) BN

(3.295)

The reﬂectance R1,n and source function ϵ+1,n regarded from the ﬁrst to the n layers as asingle layer are92

R1,n = Rn + Tn (1 − R1,n−1 Rn )−1 R1,n−1 Tn(3.296)
−1+ + T (1 − R− + ϵ+
ϵ+
=
ϵ
R
)
R
ϵ
n
1,n−1
n
1,n−1
n
n
1,n−1
1,N

(3.297)+
It can be solved by n = 2, . . . , N , starting from R1,1 = R1 , ϵ+
1,1 = ϵ1 .
With these, downward ﬂux at the boundary between layers n and n + 1, the downward andupward ﬂux are came back to a problem between two layers of combined layer, the combinationsof layers 1 − n and n + 1 − N .



−1
+
−
u+
=
(1
−
R
R
)
ϵ
+
R
ϵ
1,n
1,n
n+1,N
1,n
n+1,N
n+1/2

(3.298)+−u−
n+1/2 = Rn+1,N un,n+1 + ϵn+1,N(3.299)can be written as above. However, the ﬂux at the top of the atmosphere isu+
1/2 = 0(3.300)−u−
1/2 = ϵ1,N

(3.301)

Finally, since this ﬂux is scaled, we rescaled and added the direct solar incidence to the ﬁnd
the radiation ﬂux.
∗+−⟨τ ⟩/µ0 F
Fn+1/2= WW u+0
n+1/2 + µ0 e

(3.302)−
Fn+1/2= WW u−n+1/2(3.303)+++Also, surface direct downward radiation ﬂux Fsf
is

∗

+
Fsf
= µ0 F0 e−⟨τ ⟩ µ0

3.6.12(3.304)

Adding in the FluxF± =X

wc (1 − Ccu ) F̄ ± +X

cwc Ccu F c±

(3.305)c

F ◦± =

X

wc F ◦±

c93(3.306)

If the radiation ﬂux FC± is found for each sub-channel in each layer, the wavelength-integrated
ﬂux is found by multiplying a weight wc correspondingly to a wavelength representative of the
sub-channel and adding. Ccu is the horizontal coverage of the cumulus cloud. It is divided into
the short wavelength range and long wavelength range and added together. In addition, thedownward ﬂux of a part of the short wavelength region (shorter than the wavelength of 0.7-0.8µm) at the surface is obtained as PAR (photosynthetically active radiation). Also, the radiation
ﬂux in the clear-sky condition is calculated (F ◦± ).
Also, in the shortwave region, we ﬁnd the downward radiation at the lower end of theatmosphere and the diﬀerence between the surface direct downward radiation ﬂux.
+
Fsf
=

X

wc (1 − Cu ) F̄N++1/2 +X

c+
Fsf,dif=X

wc Ccu FNc +1/2

(3.307)

c


 X


+
c
wc (1 − Ccu ) F̄N++1/2 − F̄sf+
wc Ccu FNc +1/2 − Fsfc(3.308)cThis calculation is done in SUBROUTINE:[DTRN31]
3.6.13Calculation of the Temperature Derivative of the Flux

To implicitly solve for surface temperature, calculate diﬀerential term of upward ﬂux with respect
to surface temperature dF ∓ /dTg . Therefore, we obtained the value for temperatures 1K higher
than Tg B̄ w (Tg + 1) and used it to redo the ﬂux calculation using the addition method, andthe diﬀerence from the original value is set to dF ∓ /dTg . This is a meaningful value only in the
longwave region (earth radiation region). This calculation is done in SUBROUTINE:[RADFLX] of
pradt.F.
3.6.14

Calculation of the Heating Rate

The heating rate of the nth layer Hn is calculated by using the radiation ﬂux obtained sofar. It is calculated separately for shortwave and longwave ranges, and ﬁnally add together(SUBROUTINE:[RDTND] in pradm.F).
+(Fn− − Fn− ) − Fn+ − Fn+1
Hn = −
gCp dp3.6.15(3.309)Flux of incidence and incident angle

In this section, SUBROUTINE:[SHTINS] in pradi.F is described.
The following parameters are determined using the orbital eccentricity e, with reference toBerger (1978).
β=a1 = −2√

1 − e21
1 3
2e + 8ea2 = −2 − 14 e2



94



(3.310)(1 + β)12 +β(3.311)(3.312)a3 = −2

1 3
8e



1
2 +β



(3.313)b1 = 2e − 14 e3(3.314)b2 = 54 e2(3.315)3b3 = 1312 e(3.316)ϵ = epsd
180 π(3.317)ϖ = vpid+180π180(3.318)Additionally,

where epsd and vpid are the angle of the obliquity, and the precession represented by thetrue longitude of the perihelion measured from the vernal equinox.
The mean longitude of the vernal equinox λ0 is computed.

λ0 = a1 sin(−ϖ) + a2 sin(−2ϖ) + a3 sin(−3ϖ)(3.319)The mean longitude of the Earth position λm at time tm is represented by using the positionof the vernal equinox λ0 .
tm −t0λm = 2π×t
+ λ0
year

(3.320)where tyear is the number of seconds in year, and the origin of the time t0 is deﬁned as thetime of the vernal (March) equinox.
The true longitude of the Earth position V at time tm is calculated as below.

(3.321)V = λm − ϖ + b1 sin (λm − ϖ) + b2 sin 2 (λm − ϖ) + b3 sin 3 (λm − ϖ)(3.322)The solar declination δ is

δ = arcsin(sin ϵ sin(V + ϖ))(3.323)The incident angle cos ζ is founded by using the latitude φ, the solar declination δ, and thehour angle at a point of longitude h.
95µ0 = cos ζ = sin φ sin δ + cos φ cos δ cosh(3.324)

Incident Flux F0 is represented as follows,
F0 = F00 r−221−er = 1+e(cosV +ϖ)(3.325)(3.326)

where F00 is the solar constant and is the ratio to the time of the distance between the sun
and the earth. The number of times when cos ζ ≥ 0 (in the daytime) in time increments (set in
NHSUB), is counted, and F0 and cos ζ are ﬁnally averaged.
It is also possible to give average annual insolation. In this case, the annual and day mean
incidence and angle of incidence are approximated as follows.

3.6.16F̄ = F00 /π(3.327)µ̄0 ≃ 0.410 + 0.590 cos2 φ(3.328)

Reading Each Parameter

In SUBROUTINE:[OPPARM2]of pradt.F, various parameters used for radiation calculation are read.
The outline of the procedure is shown below.
1. Read the numbers of bands, the radiances representative of upward and downward radiation, the grids of the pressure log(p), grid of the temperatures, the optical ﬂag, and
CFCs.
2. Read the band boundaries and the information of the pressure grid and temperature grid.
3. First, the optical property ﬂag, the number of channels, the weights for channels and thenumber of the molecules including in a waveband are read. The optical property ﬂag is
modiﬁed in order to distinguish the PAR, 0.67, 0.5, and 10.5 µm. Additionally, molecule
ID and k-width are read for each molecule. Finally, the absorption coeﬃcient for the
self-broadening and CFC are also done only when the optical property ﬂag in the band ispositive. The k-width and the absorption coeﬃcient for the self-broadening are arrangedin the order of the grid of the temperatures, the grids of the pressures, and the channels.
Step 3 is performed for each wavelength band.
4. First, the number of particles is read. Next, the numbers of the modes and the mode
radius (or relative humidity) are read for each particle. Using the mode radius (or relative humidity), the following parameter required to interpolate the calculated values iscalculated for each mode number.

1(ref )(ref )
Rn+1 − Rn

96(3.329)

5. Read the band boundaries again.
6. Read the Plank function coeﬃcient, solar insolation, surface properties (not output),Rayleigh scattering coeﬃcient, Rayleigh scattering phase function. The moment for particle scattering phase function is read in the order of the particle and the optical number
and read up to the second moment. Step 6 is performed for each wavelength band.
3.6.17

Other Notes

The calculation of the radiation is usually not done at every step. Thus, the radiation ﬂux issaved, and it is used if the time is not used for radiation calculation. As for the shortwave
radiation, using the percentage of time (time is µ0 > 0) between the next calculation time fand the solar incidence angle factor averaged over the daylight hours µ̄0 , seek the Flux F̄ ,

F =f

µ0
F̄
µ̄0

97

(3.330)3.7

Turbulence Scheme

A turbulence scheme represents the eﬀect of subgrid-scale turbulence on the grid-mean prognostic variables. It calculates the vertical diﬀusion of momentum, heat, water and other tracers.
The Mellor-Yamada-Nakanishi-Niino scheme (MYNN scheme; Nakanishi 2001; Nakanishi andNiino 2004) has been used as a turbulence scheme in MIROC since its version 5, which is an
improved version of the Mellor-Yamada scheme (Mellor 1973; Mellor and Yamada 1974; Mellor
and Yamada 1982). Its closure level is 2.5. Level 3 is available, however it was not adopted as astandard option, since we could not gain large beneﬁts despite its much greater computationalcosts.
In the MYNN scheme, liquid water potential temperature θl and total water qw are used as
key variables, which are deﬁned as   Rd
Lv + LfLv
ps Cp
θl ≡ T −
ql −
qi
,
Cp
Cp
p

(3.331)

qw ≡ qv + ql + qi ,(3.332)where T and p are temperature and pressure; qv , ql and qi are speciﬁc humidity, liquid watercontent, and ice water content respectively; Cp and Rd are speciﬁc heat at constant pressure andgas constant of dry air respectively; Lv and Lf are latent heat of vaporization and fusion perunit mass respectively. ps is 1000hP a. These variables conserve for the phase change of water.
In the level 2.5, the scheme predicts the time evolution of twice turbulent kinetic energy as
a prognostic variable, which is deﬁned by
q 2 ≡ ⟨u2 + v 2 + w2 ⟩(3.333)where u, v, and w are velocities in zonal, meridional and vertical directions respectively.
In this chapter, uppercase letters represent grid-mean variables and lowercase counterparts thedeviation from the grid-mean. ⟨ ⟩ denotes an ensemble mean. In the Level 3, ⟨θl 2 ⟩, ⟨qw 2 ⟩, ⟨θl qw ⟩are also predicted, but we skip the details here.
The outline of the computational procedures is given as follows along with the names of thesubroutines. All the subroutines listed here are written in a Fortran source code of pvdfm.F.
1. Calculation of friction velocity and the Obukhov length
2. Calculation of buoyancy coeﬃcients [SUBROUTINE:VDFCND]
3. Calculation of stability functions of the Level 2 [SUBROUTINE:VDFLEV2]
4. Calculation of planetary boundary layer depth [SUBROUTINE:PBLHGT]
5. Calculation of master length scale [SUBROUTINE:VDFMLS]
6. Calculation of diﬀusion coeﬃcients, vertical ﬂuxes and their derivatives [SUBROUTINE:VDFLEV3]
7. Calculation of production and dissipation terms of twice turbulent kinetic energy [SUBROUTINE:VDFLEV3]
8. Calculation of tendencies of prognostic variables with implicit scheme
3.7.1

Surface Layer

The friction velocity u∗ and the Obukhov length LM are given as

98


1
4
u∗ = ⟨uw⟩g 2 + ⟨vw⟩g 2 ,LM = −(3.334)

Θv,g u∗ 3,kg⟨wθv ⟩g(3.335)where the subscript g indicates the values near the surface Θv and θv denote virtual potential
temperature, k the von Kármán constant, and g the acceleration of gravity. The values of thelowest model layer is used for Θv,g .
3.7.2

Calculation of the Buoyancy Coeﬃcients

The buoyancy-production term in the prognostic equation of the twice turbulent kinetic energy
contains ⟨wθv ⟩. Following Mellor and Yamada (1982), we assume the probability distribution ofθl and qw in a given grid and rewrite this term as⟨wθv ⟩ = βθ ⟨wθl ⟩ + βq ⟨wqw ⟩.

(3.336)

However, note that unlike Mellor and Yamada (1982) and Nakanishi and Niino (2004), theprobability distribution assumed here is not Gaussian. It is triangular documented in the PDFbased prognostic large-scale condensation scheme (Watanabe et al. 2009). In this case, thebuoyancy coeﬃcients, βθ and βq are written asβθ = 1 + ϵQw − (1 + ϵ)Ql − Qi − R̃abc,(3.337)βq = ϵΘ + R̃ac,(3.338)where ϵ = Rv /Rd − 1. Rv is the gas constant for water vapor, anda=

Lv ∂Qs1+
Cp ∂T T =Tl,(3.339)T ∂Qs,
Θ ∂T T =Tl(3.340)Θ Lv[1 + ϵQw − (1 + ϵ)Ql − Qi ] − (1 + ϵ)Θ,T Cp

(3.341)

b=

c=

!−1



Ql
R̃ = R 1 − a [Qw − Qs (Tl )]
2σs


−

Ql 2
,
4σs 2

σs 2 = ⟨qw 2 ⟩ − 2b⟨θl qw ⟩ + b2 ⟨θl 2 ⟩,(3.342)(3.343)where R and Ql are cloud amount and liquid water computed from the probability distribution in the grids, respectively, and Qs is saturation water vapor.
99

3.7.3

Stability Functions for the Level 2

It is known that the Mellor-Yamada Level 2.5 scheme fails to capture the behavior of growing
turbulence realistically (Helfand and Labraga 1988). Thus, the MYNN scheme ﬁrst calculates
the twice turbulent kinetic energy of the Level2 q2 2 , and then make a correction to the diﬀusion
when q < q2 , i.e., the turbulence is in a growing phase. The stability functions of the level 2,SH2 and SM 2 , required for the calculation of q2 , are represented by
Rfc − Rf,
1 − Rf(3.344)

Rf 1 − Rf
SH2 ,Rf 2 − Rf(3.345)

SH2 = SHC

SM 2 = SM C

where Rf is the ﬂux Richardson number and calculated ashi
Rf = Ri1 Ri + Ri2 − (Ri2 − Ri3 Ri + Ri2 2 )1/2 .

(3.346)

Here, Ri is the gradient Richardson number represented by
g
Ri =
Θ


 , "


 #
∂Θl
∂Qw
∂U 2
∂V 2
βθ
+ βq
+
.
∂z
∂z
∂z
∂z

(3.347)The other symbols indicate quantities independent of the environmental ﬁeld, which are
given as follows.
SHC = 3A2 (γ1 + γ2 ),(3.348)

A1 F1
,
A2 F2

(3.349)γ1,γ1 + γ2(3.350)γ 1 − C1,
F1

(3.351)

γ1,F2(3.352)1,
2SM c(3.353)Ri2 = Rf 1 SM C ,(3.354)

SM C =

Rfc =Rf 1 = B1

R f 2 = B1

Ri1 =

100

Ri3 = 4Rf 2 SM C − 2Ri2 ,(3.355)where

A1 = B 1

1 − 3γ1,6(3.356)

A2 = A1γ 1 − C1,γ1 P r(3.357)

C1 = γ 1 −1
1,(3.358)3A1 B1 3

F1 = B1 (γ1 − C1 ) + 2A1 (3 − 2C2 ) + 3A2 (1 − C2 )(1 − C5 ),(3.359)F2 = B1 (γ1 + γ2 ) − 3A1 (1 − C2 ),(3.360)B2
2A1(1 − C3 ) +(3 − 2C2 ) ,B1
B1(3.361)γ2 =and(P r, γ1 , B1 , B2 , C2 , C3 , C4 , C5 ) = (0.74, 0.235, 24.0, 15.0, 0.7, 0.323, 0.0, 0.2).
3.7.4
3.7.4.1(3.362)Master Length Scale
Original Formulation by Nakanishi (2001)

Nakanishi (2001) proposed the following formulation for the master length scale L.
1
1
1
1
+
+
,
=
L
LS
LT
LB

(3.363)where LS , LT , LB represent length scales in the surface layer, convective boundary layer,and stably stratiﬁed layer respectively. These length scales are formulated asζ≥1 kz/3.7kz/(2.7 + ζ)
0≤ζ<1
LS =

kz(1 − α4 ζ)0.2
ζ < 0,(3.364)Z ∞qz dz
LT = α1 Z0 ∞
,
q dz
0

101

(3.365) α
 2 q/NLB =α2 + α3 (qc /LT N )1/2 q/N∞∂Θv /∂z > 0 and∂Θv /∂z > 0 and∂Θv /∂z ≤ 0,

ζ≥0
ζ<0

(3.366)

where ζ ≡ z/LM is a height normalized by the Monin-Obukhov length LM , N ≡ [(g/Θ)(∂Θv /∂z)]1/2
is the Brunt-Väisälä frequency and qc ≡ [(g/Θ)⟨wθv ⟩g LT ]1/3 is a velocity scale in the convective
boundary layer.
3.7.4.2

Modiﬁcations in MIROC

The above formulation works well when the domain of the model is limited to the planetary
boundary layer (PBL) and its surrounding area. However, if the the upper troposphere isincluded, the formulation gives inappropriate behaviors depending on the conditions: e.g. LT ,the length scale of the convective boundary layer, is used in the free atmosphere, and theturbulent kinetic energy in the free atmosphere is taken into account in the calculation of LT .
In order to avoid such misbehaviors, the top height of the convective boundary layer HP BL1/2
is estimated in MIROC and we consider that the region below h = (FH HP BL )2 + H02
is theone where the PBL-derived turbulence is dominant. Here, we adopted FH = 1.5 and H0 = 500m.
Below the altitude h, equation (3.363) is used as the master length scale, but the verticalrange of the integration in LT is modiﬁed as
Z h
qz dz
,
LT = α1 Z0 hq dz

(3.367)0

and then the master length scale above h is represented as
1
1
1
1
=
+
+
L
LS
LA Lmax(3.368)where LA = α5 q/N is a length scale of air parcel vertically transported by turbulence in astably stratiﬁed layer. α5 represents the eﬀect of dissipation set to 0.53. Lmax = 500m gives theupper limit of L.
3.7.4.3

Estimation of the Top Height of the Convective Boundary Layer

Based on Holtslag and Boville (1993), HP BL is estimated using the bulk Richardson number
RiB given as

RiB =[g/Θv (z1 )][Θv (zk ) − Θv,g ](zk − zg ),[U (zk ) − U (z1 )]2 + [V (zk ) − V (z1 )]2 + Fu u∗ 2

(3.369)where zk is the altitude of a k-th model layer from the bottom at full level, z1 the altitudeof the lowest layer at full level, zg the altitude of the surface. Fu is a dimensionless tuning
parameter, andΘv,g = Θv (z1 ) + Fb
102⟨wθv ⟩g,wm(3.370)wm = u∗ /ϕm ,

ϕm =

zs
1 − 15LM(3.371)

− 1
3,(3.372)where zs is the altitude of the surface layer assumed to be 0.1HP BL . Fb is a dimensionlesstuning parameter.
RiB is successively calculated from k = 2 upward, and then if RiB exceeds 0.5 for the ﬁrst
time, it is linearly interpolated between this level and the level immediately below it. The heightsatiﬀying RiB = 0.5 is used as HP BL . Since HP BL is necessary for the calculation of zs , we ﬁrst
calculate zs using a temporary value of HP BL = z1 − zg , from which we calculate the ﬁrst guessof HP BL . Then we use this value for the recalculation of zs , and then it is used for the ﬁnal
estimate of HP BL .
3.7.5
3.7.5.1

Calculation of Diﬀusion Coeﬃcients
Twice Turbulent Kinetic Energy of Level 2

The twice turbulent kinetic energy of the level 2, q2 2 , is calculated from the following equation, which neglects the time derivative, advection and diﬀusion terms in the prognostic equationof the twice turbulent kinetic energy.
Ps + Pb − ε = 0,(3.373)where Ps and Pb denote the production terms by shear and buoyancy respectively. ε is thedissipation term. Ps and Pb are written as
Ps = −⟨wu⟩∂U∂V− ⟨wv⟩,∂z∂z(3.374)g⟨wθv ⟩,
Θrespectively. In the level 2 of the MYNN scheme, these are represented asPb ="Ps = LqSM 2

∂U
∂z

2



2 #
,

(3.376)



∂Θl
g
∂Qw
βθ
,
Pb = LqSH2
+ βq
Θ
∂z
∂z

(3.377)

ε=

+

∂V
∂z

(3.375)q3.
B1 L(3.378)From (3.373), (3.376), (3.377), and (3.378), q22 is calculated by(
q2 2 = B1 L2

SM 2

"

∂U
∂z

2


+

∂V
∂z

2 #

103

g
+ SH2
Θ


)
∂Θl
∂Qw
βθ
.
+ βq
∂z
∂z

(3.379)

3.7.5.2

Stability Functions of the Level 2.5

When q < q2 , i.e., the turbulence is in a growing phase, the stability functions of the Level
2.5 for momentum and heat, SM and SH respectively, are calculated using α = q/q2 introduced
by Helfand and Labraga (1998) as

SM = αSM 2 ,(3.380)SH = αSH2 .

(3.381)

When q ≥ q2 , SM and SH are calculated as

SM = A1

E3 − 3C1 E4,
E2 E4 + E5 E3(3.382)

S H = A2

E2 + 3C1 E5,
E2 E4 + E5 E3(3.383)where
E1 = 1 − 3A2 B2 (1 − C3 )GH ,(3.384)E2 = 1 − 9A1 A2 (1 − C2 )GH ,(3.385)E3 = E1 + 9A2 2 (1 − C2 )(1 − C5 )GH ,(3.386)

E4 = E1 − 12A1 A2 (1 − C2 )GH ,(3.387)E5 = 6A1 2 GM ,(3.388)L2
GM = 2q"L2 g
GH = − 2
q Θ

∂U
∂z

2


+

∂V
∂z

2 #
,

(3.389)



∂Θl
∂Qw
βθ
.
+ βq
∂z
∂z

(3.390)

The above formulas appear to be diﬀerent from those in Nakanishi (2001), but are equivalentand can be computed with a smaller computational cost.
1043.7.5.3

Calculation of Diﬀusion Coeﬃcients

The diﬀusion coeﬃcients for momentum, twice turbulent kinetic energy, heat and water are
represented by

KM = LqSM ,(3.391)

Kq = 3LqSM ,(3.392)

KH = LqSH ,(3.393)

Kw = LqSH ,(3.394)

respectively.
3.7.5.4

Calculation of Fluxes

The vertical ﬂuxes for U , V , q 2 , Cp T and Qw at half levels are calculated as
Fu,k−1/2 = −ρk−1/2 KM,k−1/2

Uk − Uk−1,∆zk−1/2(3.395)Fv,k−1/2 = −ρk−1/2 KM,k−1/2

Vk − Vk−1,∆zk−1/2(3.396)q 2 k − q 2 k−1,∆zk−1/2(3.397)Fq,k−1/2 = −ρk−1/2 Kq,k−1/2

FT,k−1/2 = −ρk−1/2 KH,k−1/2 Cp Πk−1/2

Fw,k−1/2 = −ρk−1/2 Kw,k−1/2

Θl,k − Θl,k−1,∆zk−1/2

Qw,k − Qw,k−1,∆zk−1/2(3.398)(3.399)respectively, where ρ denotes density and Π the Exner function. In order to perform time
integration with an implicit scheme, the derivative of each of the vertical ﬂuxes is also calculated
as∂Fu,k−1/2∂Fv,k−1/2∂Fu,k−1/2∂Fv,k−1/2
1==−=−= ρk−1/2 KM,k−1/2
,
∂Uk−1
∂Vk−1
∂Uk
∂Vk
∆zk−1/2

105(3.400)∂Fq,k−1/2∂Fq,k−1/2
1,=−= ρk−1/2 Kq,k−1/2∂q 2 k−1∂q 2 k∆zk−1/2(3.401)∂FT,k−1/2
Πk−1/21= ρk−1/2 KH,k−1/2 Cp,∂Tk−1
Πk−1 ∆zk−1/2(3.402)∂FT,k−1/2
Πk−1/2
1= −ρk−1/2 KH,k−1/2 Cp,∂TkΠk ∆zk−1/2(3.403)∂Fw,k−1/2∂Fw,k−1/2
1=−= ρk−1/2 Kw,k−1/2,∂Qw,k−1∂Qw,k∆zk−1/2(3.404)

where ∆zk−1/2 = zk − zk−1 . The ﬂuxes for other tracers are also calculated in the same way
using Kw .
3.7.6
3.7.6.1

Calculation of Turbulent Variables
Calculation of Twice Turbulent Kinetic Energy

The prognostic equation for q 2 is expressed as
dq 2
1 ∂Fq
=−
+ 2 (Ps + Pb − ε) .
dt
ρ ∂z

(3.405)In the Level 2.5, Ps , Pb , ε are written as
"

2



∂V
∂z

2 #

Ps = LqSM

∂U
∂z

,

(3.406)g
Pb = LqSH
Θ



∂Θl
∂Qw
βθ
+ βq
,
∂z
∂z

(3.407)

ε=

+

q3
.
B1 L

(3.408)Advection terms are calculated using the tracer transport routines in the dynamical core.
The turbulence scheme calculates the time evolution by the diﬀusion, production and dissipation
terms with an implicit scheme.

1063.7.6.2

Diagnosis of Variance and Covariance

The prognostic equations for ⟨θl 2 ⟩, ⟨qw 2 ⟩, ⟨θl qw ⟩ are expressed asd θl 2∂∂Θl=−wθl2 − 2 ⟨wθl ⟩− 2εθl ,dt∂z∂z(3.409)d qw 2
∂
∂Qw
2
=−
wqw
− 2 ⟨wqw ⟩− 2εqw ,dt∂z∂z(3.410)d ⟨θl qw ⟩
∂
∂Θl
∂Qw
= − ⟨wθl qw ⟩ − ⟨wqw ⟩− ⟨wθl ⟩− 2εθq .
(3.411)dt∂z∂z∂z
In the Level 2.5, the time derivative, advection, and diﬀusion terms in these equations are
neglected assuming the following local balances.
− ⟨wθl ⟩∂Θl− εθl = 0,∂z(3.412)− ⟨wqw ⟩∂Qw− εqw = 0,∂z(3.413)∂Θl∂Qw− ⟨wθl ⟩− 2εθq = 0.
(3.414)∂z∂z
In the Level 2.5 of the MYNN scheme, − ⟨wθl ⟩, − ⟨wqw ⟩, εθl , εqw , εθq are represented as− ⟨wqw ⟩− ⟨wθl ⟩ = LqSH

∂Θl
,
∂z

(3.415)− ⟨wqw ⟩ = LqSH

∂Qw
,
∂z

(3.416)

εθl =

q
θ2 ,
B2 L l(3.417)

εqw =

q
q2 ,
B2 L w(3.418)q⟨θl qw ⟩ .
B2 L(3.419)εθq =from (3.412)-(3.419), ⟨θl 2 ⟩, ⟨qw 2 ⟩, ⟨θl qw ⟩ can be diagnosed as⟨θl ⟩ = B2 L SH
2

2



∂Θl
∂z

2
,

(3.420)

2⟨qw ⟩ = B2 L SH

∂Qw
∂z

,

(3.421)⟨θl qw ⟩ = B2 L2 SH∂Θl ∂Qw.
∂z ∂z(3.422)2

2107

3.7.6.3

Treatment in the Bottom Layer

Since the lowest model layer corresponds to the surface layer where values of physical variables rapidly change in the vertical direction, the following Monin-Obukhov similarity theory isused to accurately evaluate the vertical gradient of the variables.
∂M
u∗
=
ϕm ,
∂z
kz

(3.423)

∂Θ
θ∗
=
ϕh ,
∂z
kz

(3.424)

∂Qv
qv∗
=
ϕh ,
∂z
kz

(3.425)where M is the horizontal wind velocity for the horizontal axis aligned to the direction of
the horizontal wind in the surface layer. ϕm and ϕh are the dimensionless gradient functions for
momentum and heat respectively. θ∗ and qv∗ are the scales of potential temperature and water
vapor in the surface layer respectively, and satisfy the following relationships.
⟨wm⟩g = −u2∗ ,(3.426)⟨wθ⟩g = −u∗ θ∗ ,(3.427)⟨wqv ⟩g = −u∗ qv∗ ,(3.428)where m is the deviation of M from the grid mean. Using M and m, the production termof the turbulence kinetic energy is written as

Ps + Pb = ⟨wm⟩

∂Mg
+ ⟨wθv ⟩.
∂zΘ

(3.429)Using (3.423), (3.426) and the deﬁnition of the Obukhov length, it is rewritten as

Ps + Pb =u3∗[ϕm (ζ1 ) − ζ1 ] ,kz1(3.430)where ζ1 is ζ at the full level of the lowest model layer.
Assuming that the eﬀect of cloud particles are negligible in the surface layer, ⟨θl 2 ⟩, ⟨qw 2 ⟩,⟨θl qw ⟩ is calculated diagnostically from (3.412)-(3.414), (3.417)-(3.419), (3.424), (3.425), (3.427),
and (3.428) as
ϕh (ζ1 )⟨wθ⟩g 2⟨θl ⟩ =u∗ kz1
2

108



q
,
B2 L

(3.431)

⟨qw 2 ⟩ =ϕh (ζ1 )⟨wqv ⟩g 2
u∗ kz1ϕh (ζ1 )⟨θl qw ⟩ =⟨wθ⟩g ⟨wqv ⟩gu∗ kz1

q
,
B2 L


q
.
B2 L

(3.432)(3.433)ϕm and ϕh are formulated following Businger et al. (1971) as

1 + β1 ζ,ζ≥0−1/4(1 − γ1 ζ), ζ<0(3.434)β2 + β1 ζ,ζ≥0−1/2
β2 (1 − γ2 ζ), ζ<0(3.435)(β1 , β2 , γ1 , γ2 ) = (4.7, 0.74, 15.0, 9.0).

(3.436)ϕm (ζ) =ϕh (ζ) =3.7.7

Time Integration with Implicit Scheme

3.7.7.1

Tendency of q 2

The prognostic equation for q 2 is discretized as22
qk,n+1− qk,n

∆t

!
turb




qk,n 21(3.437),F− Fq,k−1/2,n+1 + 2 Ps,k,n + Pb,k,n −q=−ρk ∆zk q,k+1/2,n+1
B1 L k,n+1

where n and n + 1 indicate the current and next time steps respectively, and ∆zk ≡ zk+1/2 −zk−1/2 . The subscript turb indicates the calculation by the turbulence scheme and the advection
term is omitted. Fq at n + 1 is computed by

Fq,k−1/2,n+1 = Fq,k−1/2,n +∂Fq,k−1/2 2∂Fq,k−1/2 2
2
2(qk,n+1 − qk,n)+(qk−1,n+1 − qk−1,n). (3.438)
2
2
∂qk
∂qk−1

With a deﬁnition of
µk =

2
2
qk,n+1− qk,n∆t!
,(3.439)turb(3.437) and (3.438) lead to

X1,k µk+1 + X2,k µk + X3,k µk−1 = Yk ,
where

109

(3.440)∂Fq,k+1/2
∆t,
2
∂qk+1

 

∂Fq,k+1/2 ∂Fq,k−1/2
qk,nX2,k = ρk ∆zk 1 + 2
∆t +
−
∆t,
B1 L
∂qk2
∂qk2
∂Fq,k−1/2
X3,k = −
∆t,
2
∂qk−1
X1,k =

Yk = Fq,k−1/2,n − Fq,k+1/2,n + 2ρk ∆zk

Ps,k,n + Pb,k,n −3qk,n(3.441)(3.442)(3.443)!

B1 L.

(3.444)(3.440) makes the following matrix equation,


X2,K
 X1,K−1

 0

 ..
 .

 0

 0
0

X3,K
0
0
0
X2,K−1 X3,K−1 00
X1,K−2 X2,K−2 X3,K−2 0
..
.
···
···
···

0
0
0

X1,3
0
0

···
···
···

X2,3 X3,3
X1,2 X2,2
0
X1,1



 

µK
YK
  µK−1   YK−1    µK−2   YK−2    ..
  ..

,
 .
= .
(3.445)

 






0
  µ3  Y3  Y2
X3,2   µ2
X2,1
µ1
Y1

0
0
0
..
.

where the subscript K denote the index for the top model layer. (3.445) is solved for µk
using the LU decomposition.
3.7.7.2

Tendencies of the Other Prognostic Variables

Letting ψ be a substitute for u, v, T , qw , the tendency of ψ is calculated byψk,n+1 − ψk,n
∆t


=−
turb


1
Fψ,k+1/2,n+1 − Fψ,k−1/2,n+1 ,ρk ∆zk

(3.446)whereFψ,k−1/2,n+1 = Fψ,k−1/2,n +∂Fψ,k−1/2∂Fψ,k−1/2(ψk,n+1 − ψk,n ) +(ψk−1,n+1 − ψk−1,n ).(3.447)
∂ψk
∂ψk−1

These equations lead to (3.445) again and computed with the LU decomposition, but µk ,X1,k , X2,k , X3,k and Yk are redeﬁned asψk,n+1 − ψk,n
µk =
,
∆t
turb
∂Fψ,k+1/2
X1,k =
∆t,
∂ψk+1


∂Fψ,k+1/2 ∂Fψ,k−1/2∆t,X2,k = ρk ∆zk +
−
∂ψk
∂ψk
∂Fψ,k−1/2
X3,k = −∆t,∂ψk−1
Yk = Fψ,k−1/2,n − Fψ,k+1/2,n .
110(3.448)(3.449)(3.450)(3.451)(3.452)

3.8

Surface Flux Scheme (Sea Surface)

Until CCSR/NIES AGCM (1997), both land surface and sea surface were treated as one of
the atmospheric physical processes, but after MIROC3 (Hasumi and Emori, 2004), land surface processes became independent as MATSIRO. However, since MIROC3 (Hasumi and Emori,2004), land surface processes have been separated into MATSIRO (Takata et al., 2003; Nitta etal., 2014). In SUBROUTINE:[SURFCE] in pgsfc.F, ENTRY:[OCNFLX] (in SUBROUTINE:[OCEAN]
of pgocn.F) is called for the sea surface, and ENTRY:[LNDFLX] (in SUBROUTINE:[MATSIRO]of matdrv.F) is called for the land surface, respectively. This chapter describes sea surface
processes, which are still treated within the framework of atmospheric physical processes in
MIROC6 (Tatebe et al., 2019)). For the land surface processes, please refer to Description of
ILS (https://github.com/integrated-land-simulator/model description).
Sea surface processes provide the boundary conditions at the lower end of the atmospherethrough the exchange of momentum, heat, and water ﬂuxes between the atmosphere and thesurface. In ENTRY:[OCNFLX], the following procedure is used to deal with sea surface processes.
1. Prepare variables for sea ice extent and no ice extent, respectively, using sea ice concentration.
2. Determine the surface boundary conditions.
3. Calculate the ﬂux balance.
4. Calculate the radiation budget at the sea surface.
5. Calculate the deposition by CHASER.
6. Solve the heat balance at the sea surface and update the skin temperature and each ﬂux
value.
No prognostic variables are used in this scheme.
Practically, precipitation ﬂux from 2 schemes are treated together.

P r = P rc + P rl

(3.453)where P r is total precipitation ﬂux, P rc is precipitation ﬂux from the cumulus convection
scheme, and P rl is precipitation ﬂux from the large scale condensation scheme, respectively.
Sea ice covered/free areas are represented by L = 1, 2. Each area is calculated then weighted
by sea ice concentration (Rice ).
In the sea ice area (L = 1), the skin temperature (Ts ) is the sea ice skin temperature (Tice ).
However, if Tice is higher than the melting point (Tmelt = 273.15[K]), then Tmelt is used.
Ts = min(Tice , Tmelt )(3.454)The sea ice bottom temperature (Tb ) is assumed to be the sea skin temperature (To(1) ).
Tb = To(1)(3.455)The amount of sea ice (Wice ) and the amount of snow on it (Wsnow ) are converted per unit
area by considering sea ice concentration (Rice ) and used in the calculation. However, a limiter(ϵ) is provided to prevent the values from becoming too small.

Rice = max(Rice,orginal , ϵ)111

(3.456)

In the ice-free area (L = 2), the skin temperature (Ts ) and sea ice bottom temperature (Tb )are assumed to be the sea temperature temperature (To(1) ).
Ts = Tb = To(1)(3.457)The evaporation eﬃciency is set to 1 for both L = 1, 2.
If sea ice concentration (Rice ) is not given (as a boundary condition or from an OGCM), it issimply diagnosed with the sea ice volume (Wice ) in ENTRY:[OCNICR] (in SUBROUTINE:[OCNICR]
of pgocn.F).
sRice = minmax(Wice , 0)
, 1.0
Wice,c(3.458)The standard gives the amount of sea ice per area as Wice,c = 300[kg/m2 ].
3.8.1

Boundary Conditions

In ENTRY[OCNBCS] (in SUBROUTINE:[OCNSUB] of pgocn.F), surface albedo and roughness arecalculated. They are calculated supposing ice-free conditions, then modiﬁed to take into accountthe eﬀects of ice and snow cover.
First, let us consider the sea surface level albedo (α(d,b) ), b = 1, 2, 3 represent the visible,
near-infrared, and infrared wavelength bands, respectively. Also, d = 1, 2 represents direct andscattered light, respectively. The albedo for the visible bands are calculated in SUBROUTINE[SEAALB] (of pgocn.F), supposing ice-free conditions. The albedo for near-infrared is set tosame as the visible one. The albedo for infrared is uniformly set to a constant value.
The grid-averaged albedo, taking into account sea ice concentration (Rice ), isα = α − Rice αice(3.459)αice is given as αice,1 = 0.5, αice,2 = 0.5, αice,3 = 0.05, respectively.
In addition, we want to consider the eﬀect of snow cover. Here, we consider the albedo
modiﬁcation by temperature. Standard threshold values for snow temperature are Tal,2 =258.15[K] and Tal,1 = 273.15[K]. The snow albedo changes linearly with temperature change
from αsnow,1 = 0.75 to αsnow,2 . Let the coeﬃcient τsnow , which is 0 ≤ τ ≤ 1.

τsnow =

Ts − Tal,1
Tal,2 − Tal,1(3.460)Update the snow albedo (αsnow ) as
αsnow = αsnow,0 + τsnow (αsnow,2 − αsnow,1 )(3.461)

Second, let us consider sea surface roughnesses. Roughnesses of for momentum, heat andvapor are calculated in SUBROUTINE:[SEAZ0F] (of pgocn.F), supposing the ice-free conditions,then modiﬁed to take into account the eﬀects of ice and snow cover.
When the sea ice exists (L = 1), roughnesses of momentum, heat and vapor (r0,M , r0,H , r0,E )is modiﬁed to take into account sea ice concentration (Rice ),112

z0,M = z0,M + (z0,ice,M − z0,M )Rice(3.462)z0,H = z0,H + (z0,ice,H − z0,H )Rice(3.463)z0,E = z0,E + (z0,ice,E − z0,E )Rice(3.464)

where, r0,ice,M , r0,ice,H , r0,ice,E is the roughness of sea ice for momentum, heat and vapor,respectively.
z0,M = z0,M + (z0,snow,M − z0,M )Rsnow(3.465)z0,H = z0,H + (z0,snow,H − z0,H )Rsnow(3.466)z0,E = z0,E + (z0,snow,E − z0,E )Rsnow(3.467)

where, r0,snow,M , r0,snow,H , r0,snow,E is the roughness of snow ice for momentum, heat andvapor, respectively.
Third, let us consider conductivity of ice. When sea ice exists (L = 1), thermal conductivity⋆ ) is obtained byof sea ice (kice
⋆
kice
=

Df,icemax(Rice /σice , ϵ)(3.468)where Df,ice is thermal diﬀusivity of sea ice, and σice is sea ice density, respectively.
The calculated thermal conductivity is modiﬁed to kice to take into account that it varies
with snow cover.
hsnow = min(max(Rsnow /σsnow ), ϵ), hsnow,max )⋆kice = kice
(1 − Rice ) +

Dice
Rice
1 + ∥Dice /Dsnow · hsnow ∥(3.469)(3.470)where hsnow is snow depth, Rsnow is snow area fraction, σsnow is snow density, hsnow,max is
maximum snow depth, and Dsnow is thermal diﬀusivity of snow, respectively.
Therefore, heat conduction ﬂux and its derivative are
G = kice (Tb − Ts )∂G
= kice∂T
Note that in the ice-free area (L = 2)
G = kocn(3.471)(3.472)(3.473)where kocn is heat ﬂux in the sea temperature layer, and kocn is heat ﬂux in the sea temperature layer, respectively.
113

3.8.1.1

Albedo for Visible

In SUBROUTINE [SEAALB] (of pgocn.F), albedo for the visible bands are calculated supposing
ice-free conditions, then modiﬁed to take into account the eﬀects of ice and snow cover.
For sea surface level albedo (αL(d) ), d = 1, 2 represents direct and scattered light, respectively.
Using the solar zenith angle at latitude ζ, albedo for direct light is presented by∗∗αL(1) = e(C3 A +C2 )A +C1(3.474)A∗ = min(max(cos(θ), 0.03459), 0.961)(3.475)

where

and C1 , C2 , C3 are constant parameters, respectively.
On the other hand, albedo for scattered light (αL(2) ) is uniformly set to a constant parameter.
αL(2) = 0.06
3.8.1.2(3.476)

Roughnesses

In SUBROUTINE:[SEAZ0F] (of pgocn.F), the roughnesses of for momentum, heat and vapor
are calculated supposing the ice-free conditions. calculated, according to Miller et al. (1992),
then modiﬁed to take into account the eﬀects of ice and snow cover.
The roughness variation of the sea surface is determined by the friction velocity (u⋆ ).

u⋆ =qCM0 (ua 2 + va 2 )(3.477)

where CM0 is a bulk coeﬃcient for momentum, and ua , va are zonal and vertical winds of the
1st layer of the atmosphere. We perform successive approximation calculation of CM0 , because
Fu , Fv , Fθ , Fq are required.
Then, roughnesses of sea surface for momentum, heat and vapor are

z0,M = z0,M0 + z0,MR +z0,MR u⋆ 2 z0,MS ν+gu⋆(3.478)z0,H = z0,H0 + z0,HR +z0,HR u⋆ 2 z0,HS ν+gu⋆(3.479)z0,E = z0,E0 + z0,ER +z0,ER u⋆ 2 z0,ES ν+gu⋆(3.480)where, ν = 1.5 × 10−5 [m2 /s] is the kinetic viscosity of the atmosphere, z0,M , z0,H and z0,E
are surface roughness for momentum, heat, and vapor, z0,M0 , z0,H0 and z0,E0 are base, and roughfactor (z0,MR , z0,MR and z0,ER ), and smooth factor (z0,MS , z0,MS and z0,ES ), respectively.

114

3.8.2

Calculation of Momentum, Heat and Water Vapor Fluxes

Treatment of sea surface ﬂux is basically the same with CCSR/NIES AGCM (1997). The
surface ﬂux scheme evaluates the physical quantity ﬂuxes between the atmospheric surfaces due
to turbulent transport in the boundary layer. The main input are horizontal wind speed (ua , va ),temperature (Ta ), and speciﬁc humidity (qa ) from the 1st layer of the atmosphere. The output
are the vertical ﬂuxes and the diﬀerential values (for obtaining implicit solutions) of momentum,
heat, and water vapor.
Surface ﬂuxes (Fu , Fv , Fθ , Fq ) are expressed using bulk coeﬃcients for momentum, head andvapor (CM , CH , CE ) as follows
Fu = −ρCM |Va |ua(3.481)

Fv = −ρCM |Va |va(3.482)Fθ = ρcp CH |Va |(θs − θa )(3.483)

FqP = ρCE |Va |(qs − qa )(3.484)

where FqP is the possible evaporation ﬂux, where Va is horizontal wind vector, and θs , θa arepotential temperature of surface and 1st layer of the atmosphere, respectively. Although there
is no description, surface ﬂuxes are calculated using wind speed relative to ocean current speed.
For example, F u=-￥rho C M | V a - V o | (u a - u o). Here, V o = ( u o, v o) represents
ocean current vector at the uppermost layer. Note that in a stand-alone AGCM, V o = (0,0) is
assumed.
Turbulent ﬂuxes at the sea surface are solved by bulk formulae as follows. Then, by solving
the surface energy balance, the ground skin temperature (Ts ) is updated, and the surface ﬂux
values with respect to those values are also updated. The solutions obtained here are temporary
values. In order to solve the energy balance by linearizing with respect to Ts , the diﬀerential
with respect to Ts of each ﬂux is calculated beforehand.
• Momentum ﬂux

τx = −ρCM |Va |ua(3.485)τy = −ρCM |Va |va(3.486)

where τx and τy are the momentum ﬂuxes (surface stress) of the zonal and meridional
directions, respectively.
• Sensible heat ﬂux

115

Hs = cp ρCHs |Va |(Ts − (Ps /Pa )κ Ta )(3.487)

where Hs is the sensible heat ﬂux from the sea surface; κ = Rair /cp and Rair are the gas
constants of air, and cp is the speciﬁc heat of air.
• Bare sea surface evaporation ﬂux

FqP = ρCE |Va | (q ∗ (Ts ) − qa )3.8.2.1(3.488)

Bulk factors

In SUBROUTINE:[BLKCOF] (of psfcl.F), the bulk factors are calculated. The bulk Richardson
number (RiB ), which is used as a benchmark for the stability between the atmospheric surfaces,isgg Ta (ps /pa )κ − T0θs (θa − θ(z0 ))/za
RiB =
=
fT
(ua /z1 )2θsu2a /z1(3.489)Here, g is the gravitational accerelation, θs is surface potential temperature, Ta is the atmospheric temperature of the 1st layer, Ts is the surface skin temperature, ps is the surfacepressure, pa is the pressure of the 1st layer, κ is the Karman constant, andfT = (θa − θ(z0 ))/(θa − θs )(3.490)The bulk coeﬃcients of CM , CH , CE are calculated according to Louis (1979) and Louis etal. (1982) However, corrections are made to take into account the diﬀerence between momentum and heat roughness. If the roughnesses for momentum, heat, and water vapor are set toz0,M , z0,H , z0,E , respectively, the results are generally z0,M > z0,H , z0,E , but the bulk coeﬃcientsgffor heat and water vapor for the ﬂuxes from the height of z0,M are also set to C
H , CE , thencorrected.

CM =


−eM

iB ]
 C0,M [1" + (bM /eM )R C0,M 1 − bM RiB 1 + dM bM C0,M

g
C
H =


−eH
]

0,H [1iB ] C" + (bH /eH )R]]

0,H 1 − bH RiB 1 + dH bH C0,H
 C

q

q

z1
z0,M |RiB |z1z0,M |RiB |

−1 #

−1 #

, RiB ≥ 0
, RiB < 0

, RiB ≥ 0
, RiB < 0g
CH = C
H fT

fE =
C


g + (bE /eE )RiB ]−eE
 C0,E [1
"


g
g

 C0,E 1 − bE RiB 1 + dE bE C0,E
116

(3.491)(3.492)(3.493)qz1
z0,M |RiB |−1 #, RiB ≥ 0
, RiB < 0(3.494)

fE fqCE = C(3.495)gg
C0M , C
0H , C0E is the bulk coeﬃcient (for ﬂuxes from z0M ) at neutral,k2gg
C0M = Ci2
0H = C0E = h 1ln zz0M(3.496)

fq = (qa − q(z0 ))/(qa − q ∗ (θ0 ))(3.497)

Correction Factor fq is ,but the method of calculation is omitted. The coeﬃcients of Louis factors are (bM , dM , eM ) =(9.4, 7.4, 2.0), (bH , dH , eH ) = (bE , dE , eE ) = (9.4, 5.3, 2.0).
is a correction factor, which is approximated from the uncorrected bulk Richardson number,
but we abbreviate the calculation here.
3.8.3

Radiation Flux Calculation

In SUBROUTINE:[RADSFC] (of pgsfc.F), the radiation ﬂux at sea surface is calculated. For the
ground surface albedo (α(d,b) ), b = 1, 2 represent the visible and near-infrared wavelength bands,
respectively. Also, d = 1, 2 are direct and scattered, respectively. For the downward shortwave
radiation (SW ↓ ) and upward shortwave radiation (SW ↑ ) incident on the earth’s surface, thedirect and scattered light together are↓↓↓↓SW ↓ = SW(1,1)+ SW(1,2)+ SW(2,1)+ SW(2,2)↓↓↓↓
SW ↑ = SW(1,1)· α(1,1) + SW(1,2)· α(1,2) + SW(2,1)· α(2,1) + SW(2,2)· α(2,2)

3.8.4(3.498)

Solving Heat Balance

In SUBROUTINE:[OCNSLV] (of pgocn.F), heat balance at the sea surface is solved. Downward
radiative ﬂuxes are not directly dependent on the condition of the sea surface, and their observed
values are simply speciﬁed to drive the model (Hasumi, 2015). Shortwave emission from the sea
surface is negligible, so the upward part of the shortwave radiative ﬂux is accounted for solely
by reﬂection of the incoming downward ﬂux. Let αS be the sea surface albedo for shortwave
radiation. The upward shortwave radiative ﬂux is represented bySW ↑ = −αS SW ↓(3.499)On the other hand, the upward longwave radiative ﬂux (LW ↑ ) has both reﬂection of theincoming ﬂux and emission from the sea surface. Let α be the sea surface albedo for longwave
radiation and ϵ be emissivity of the sea surface relative to the black body radiation, respectively.
The upward shortwave radiative ﬂux is represented by
LW ↑ = −αLW ↓ + ϵσTs4117(3.500)where σ is the Stefan-Boltzmann constant and Ts is skin temperature, respectively . If sea
ice exists (L = 1), snow or sea ice temperature is considered by fractions. When radiativeequilibrium is assumed, emissivity becomes identical to co-albedo:
ϵ=1−α(3.501)

The net surface ﬂux (F ∗ ) is presented byF ∗ = H + (1 − αLk )σTs4 + αLk LW ↓ − LW ↓ + SW ↑ − SW ↓(3.502)where H is sensible heat ﬂux.
With the surface heat ﬂux calculated in SUBROUTINE:[SFCFLX] (of psfcm.F) (G), heat ﬂux
into the sea surface (G∗ ) is presented asG∗ = G − F ∗(3.503)Note that G∗ is downward positive.
The temperature derivative term of G∗ is
∂G∗
∂G
∂H
∂R
=
+
+
∂Ts
∂Ts ∂Ts ∂Ts(3.504)

When the sea ice exists (L = 1), the surface ﬂux Gice is considered with the sublimation ﬂux(ls E).
Gice = G∗ − ls E

(3.505)The temperature derivative term of Gice is
∂Gice
∂G∗
∂E
=
+ ls
∂Ts
∂Ts
∂Ts

(3.506)ice −1
We can update the skin temperature with sea ice concentration and ∆Ts = Gice ( ∂G∂Ts )Ts = Ts + Rice ∆Ts(3.507)

Then, the sensible and latent heat ﬂux on the sea ice (Eice , Hice ) is updated.

Eice = E +

∂E
∆Ts
∂Ts

(3.508)Hice = H +

∂H
∆Ts
∂Ts

(3.509)When the sea ice does not existed (L = 2), otherwise, the surface heat ﬂux (Gf ree ) is
calculated by addition of evaporation ﬂux lc E and the net ﬂux F ∗ .
118

Gf ree = F ∗ + lc E

(3.510)

Finally each ﬂux is updated. For sensible heat ﬂux (H), the temperature change on the sea
ice is considered.

H = H + Rice Hice

(3.511)

Then, the heat used for the temperature change (F ) is saved.

F = Rice Hice

(3.512)For upward longwave radiative ﬂux (LW ↑ ), temperature change on the sea ice (∆Ts ) is
considered.
LW ↑ = LW ↑ + 4

σ
Rice ∆Ts
Ts(3.513)For the surface heat ﬂux (G), sea ice existence is considered.
G = (1 − Rice )Gf ree + Rice Gice

(3.514)For latent heat ﬂux E, sea ice existence is considered.
E = (1 − Rice )E + Rice Eice(3.515)

Then, each term above are saved as freshwater ﬂuxes (Wf ree , Wice ) of ice covered and freeareas.
Wf ree = (1 − Rice )E(3.516)

Wice = Rice Eice(3.517)

119

3.9

Gravity Waves

Eﬀects of sub-grid scale gravity waves on the grid scale ﬂows are represented by two kinds ofgravity wave parameterizations. The orographic gravity wave drag scheme represents decelerations of the grid scale ﬂow due to momentum deposition of orographic gravity waves, which is
described below. In contrast, the non-orographic gravity wave scheme may either accelerate ordecelerate the grid scale ﬂow. At present this document does not contain a detailed descriptionof the non-orographic gravity wave scheme, and readers may ﬁnd if in Watanabe (2008) and
Tatebe et al. (2019).
3.9.1

Overview of a Gravity Wave Drag Parameterization

The orographic gravity wave drag scheme represents the upward momentum ﬂux of the gravity
waves induced by sub-grid scale topography and calculates the horizontal wind deceleration
associated with its convergence (McFarlane, 1987). The main input data are eastward wind (u),northward wind (v), and temperature (T ), and the output data are the tendency of eastward
wind and northward wind, ∂u/∂t, ∂v/∂t.
The outline of the calculation procedure is as follows.
1. The momentum ﬂux at the ground surface is calculated from the variance of sub-grid scale
orography, the horizontal wind speed at the lowest level, and the static stability.
2. The upward propagation of gravity waves is considered. If the momentum ﬂux exceeds
its critical value, which is determined by the critical Froude number, then wave breaking
occurs and the momentum is limited by the critical value.
3. The tendency of horizontal wind is obtained by calculating the vertical convergence of
momentum ﬂux in each layer.
3.9.2

Relationship between Local FroudeNnumber and Momentum Flux

Considering the vertical ﬂux of horizontal momentum due to orographic gravity waves, thediﬀerence between the ﬂux (τ ) and the local Froude number (FL = N H/U ) at a certain altitude
is

FL =

τN
Ef ρU 3

1/2,

(3.518)This relationship holds for the following cases where N = g/θ∂θ/∂z is the Brant-Väisälä
frequency, ρ is the density of the atmosphere, U is the wind speed, and Ef is the proportional
constant corresponding to the horizontal scale of the rippling of the surface. From now on,

τ=

Ef FL2 ρU 3
N

(3.519)The local Froude number (FL ) cannot exceed the critical Froude number (Fc ). If the local
Froude number calculated from (3.518) exceeds the critical Froude number Fc , the gravity wave
becomes supersaturated and the ﬂux decreases to the momentum ﬂux corresponding to thecritical Froude number.

120

3.9.3

Momentum Fluxes at the Surface

The magnitude of the vertical ﬂux of horizontal momentum due to gravity waves excited at theearth’s surface, τ1/2 , is calculated by substituting the local Froude number (FL )1/2 = N1 h/U1
into (3.519),
τ1/2 = Ef h2 ρ1 N1 U1 ,(3.520)where U1 = |v1 | = (u21 + v12 )1/2 is the surface wind speed, N1 , ρ1 are estimated to be thestability and density of the atmosphere near the earth’s surface, respectively. h is an indicatorof the change in the sub-grid orography and is assumed to be equal to the standard deviationof the surface height (ZSD ).
Here, when the local Froude number ((FL )1/2 = N1 ZSD /U1 ) exceeds the critical Froudenumber (Fc ), the momentum ﬂux is suppressed to the value obtained by substituting (3.519) for
Fc . In other words,
Ef Fc2 ρ1 U13
2
τ1/2 = min Ef ZSD ρ1 N1 U1 ,
N1
3.9.4

(3.521)

Momentum Fluxes in the Upper Levels

Suppose that the momentum ﬂux τk−1/2 is computed at level k − 1/2. When no saturation
occurs, τk+1/2 is equal to τk−1/2 . If the momentum ﬂux (τk−1/2 ) exceeds the momentum ﬂux
calculated from the critical Froude number at the k + 1/2 level, wave breaking occurs in the klayer and the momentum ﬂux decreases to the critical ﬂux.

τk+1/2 = min τk−1/2 ,3
Ef Fc2 ρk+1/2 Uk+1/2

Nk+1/2!
,(3.522)

Note that Uk+1/2 is the magnitude of the horizontal wind speed projected on to the directionof the lowest level of the horizontal wind,

Uk+1/2 =3.9.5

vk+1/2 · v1
|v1 |(3.523)The Magnitude of the Time Variation of Horizontal Wind due to Momentum
Convergence

The tendency of the projected component of the horizontal wind, Uk , is ,
∂U
1 ∂τ
∂τ
=−
=g .
∂t
ρ ∂z
∂p

(3.524)τk+1/2 − τ k − 1/2
∂Uk
=g
.
∂t
∆p

(3.525)That isUsing this, tendency for the eastward and northward winds are calculated as follows
121

∂uk
∂Uk u1
=
∂t
∂t U1
∂vk
∂Uk v1
=
∂t
∂t U1
3.9.6

(3.526)(3.527)

Other Notes

1. It is assumed that no gravity waves are excited at the ground surface when the wind speed
is small (U1 ≤ vmin ) or when the undulations at the surface are small (ZSD ≤ (ZSD )min ).

122

4

Miscellaneous

4.1

Coupler

4.1.1

Fluxes to Atmospheric Models

4.1.1.1

Fluxes between Atmosphere and Ocean

Fluxes from the sea surface to the atmosphere (F LXO) are calculated on the sea surface
grid of the atmospheric model.
Boundary conditions such as sea surface temperature, sea ice concentration, sea ice thickness,snow depth over sea ice, sea ice internal temperature, and ocean surface current velocity are
obtained from the ocean model through an exchanger (sea ice surface temperature is determinedfrom the sea ice internal temperature, sea ice thickness, and atmospheric conditions over the sea
ice. The sea ice velocity is not currently used for ﬂux calculations in MIROC). The atmosphericboundary conditions such as wind speed, temperature, and speciﬁc humidity at sea are converted
from the atmospheric grid to the sea surface grid using linear or cubic spline completion. The
ﬂuxes from the sea surface are calculated separately for seawater and sea ice, averaged by area
weight, and passed to the atmosphere. When using a sea ice model categorized by sea ice
thickness, it may be necessary to calculate ﬂuxes for each sea ice thickness category, but thecurrent model speciﬁcation calculates ﬂuxes for the average sea ice thickness. The conversionof ﬂuxes and boundary conditions between the atmosphere and ocean by the exchanger will be
described in detail later.
4.1.1.2

Fluxes between Atmospheric land Surface

Fluxes from the land surface to the atmosphere (F LXL) are calculated on a land surface
grid. A land surface grid consists of multiple soil covers and lakes. Freezing and thawing oflakes and snow cover are considered by a vertical 1D ice model (0-layer model). If the area ofthe land surface grid is SL, the area occupied by the lake and each soil cover is respectively
SLlake = SL × LKF RC × F LN D

(4.1)

SLgrd
k = SL × GRF RCk × (1 − LKF RC) × F LN D

(4.2)

where LKF RC is the percentage of lakes on land, k is the type of soil cover, and GRF RCk is
the percentage of soil cover k on land excluding lakes. Fluxes from the land surface are calculated
separately over each of these soil covers and lakes, averaged by area weight, and passed to theatmosphere.

F LXL = LKF RC × F LXLlake+ (1 − LKF RC) ×km
X(GRF RCk × F LXLgrd
k )

(4.3)

k=1

where F LXLlake is the ﬂux at the lake surface, F LXLgrd
is the ﬂux at soil cover k, and kmkis the number of soil cover types.

1234.1.1.3

Total Flux to the Atmosphere

Since the river is treated without area property in the model, the ﬂux to the atmosphere(F LXA) can be obtained as a weighted average of the sea-land distribution of the ﬂuxes on theland grid (F LXL) and at the sea grid (F LXO) as follows
XX
1
land
× F LXLij )F LXA =
×[
(SLij × F LN Dij
SA
jldiv ildiv
j=1 i=1+

jodiv
X iodiv
X

oc(SOij × (1 − F LN Dij) × F LXOij )](4.4)

j=1 i=1.
where, (ildiv,jldiv) is the number of east-west and north-south divisions of the land surface
grid. Fluxes computed in the atmospheric model, such as precipitation, are also included in
F LXL and F LXO. In the case of such ﬂuxes, all the ﬂuxes in the partitioned land and sea
surface grids have the same value as the corresponding grid.
4.1.2
4.1.2.1

Fluxes between Land Surface Model and River Model
Fluxes between Land Surfaces and the River

In the current speciﬁcation of the model, the ﬂuxes of water between river and land surfaces
deal only with the inﬂow of water from the river to the lake (RU N IN ), the outﬂow from thelake to the river (RU N OF F ), the inﬂow of water to the land surface at the inland vanishingpoint (RU N IN ), and the outﬂow of water overﬂowing the soil to the river (RU N OF F ). Here,the inland vanishing point indicates the point where the endpoint of the river disappears, suchas in deserts. The water balance in the river model is divided into ice and water. The ice inthe river model corresponds to a pseudo-glacier. Here, the phase change in the river model is
not considered to guarantee the conservation of melting heat. In addition, the ﬂow rate of ariver is deﬁned as the amount of water present in the river grid divided by its area. Water andice are transported downstream in the river model according to the river channel network data.
In the river model in MIROC6, the river discharge at the inland vanishing point of the river isscattered over the global ocean to obtain the water balance.
4.1.2.2

Water Runoﬀ from Land Surface

When each soil cover in the land surface grid can no longer hold water or snow and ice, wateror ice is passed from each soil model to the river model through the coupler.

grd
RU N OF Fall
= (1 − LKF RC) ×km
X(GF LRCk × RU N OF Fkgrd )

(4.5)

k=1

The details of the runoﬀ from each soil cover can be found in the documentation of the land
surface model MATSIRO. In the lake model, when the lake level or snow/ice thickness (H)exceeds a constant value (Hc ), the water ﬂows out to the river at a time constant τh

RU N OF F lake = LKF RC ×

(H − Hc )
,
τh

124

(H > Hc )(4.6)

RU N OF F lake = 0,

(H < Hc )(4.7)

The average runoﬀ from the land surface is as follows.
grd
land
RU N OF Fall
= RU N OF F lake + RU N OF Fall

(4.8)

When considering the average runoﬀ volume of the land surface grid, it is necessary tomultiply the above equation by the percentage of land surface F LN D. In the river model,
land is converted to the river grid with the weight of sea-land distribution, and the
RU N OF Fallrunoﬀ amount RU N OF F riv is used for calculation.
4.1.2.3

Runin of Water from a River to a Lake

When a lake exists in the middle of a river channel, water ﬂows into the lake according to
the river ﬂow rate. In order to calculate the amount of water ﬂowing into the lake, the river
ﬂow GDRIV in the river grid is converted to the river ﬂow GDRIV L in the land surface gridthrough the coupler. Here, GDRIV L is the amount normalized by the area of the land surface
grid. In the land surface grid, the river inﬂow to the lake, RU N IN N , is deﬁned by the river
ﬂow (GDRIV L) and the time constant τ as follows
RU N IN N lake = GDRIV L/τ

(4.9)

.
Since the current speciﬁcation only considers inﬂow from rivers to lakes, except at the inland
vanishing point, the average inﬂow at the land surface is
RU N IN N land = RU N IN N lake × LKF RC(4.10).
When there are multiple river grids corresponding to a land surface grid, if the river waterinﬂow to the land surface averaged over the land surface grid is returned to the river grid using
only the area weights as in RU N OF F , it is possible that more water will ﬂow out of the riverthan exists in the river grid. Therefore, we convert the ratio of discharge to river ﬂow from theland surface grid to the river grid, and estimate the river discharge (inﬂow to the land surface)in each river grid. The runoﬀ ratio of the river ﬂow to the land surface grid is
RIN N land = RU N IN N land /GDRIV L(4.11).
If the discharge rate converted to the river grid is RIN N riv , the discharge (inﬂow to theland surface) in the river grid is
RU N IN N riv = RIN N riv × GDRIV.
125(4.12)

4.1.3Fluxes to the Ocean Model

4.1.3.1

Boundary Conditions for the Ocean on a Sea Level Grid

As mentioned above, the ﬂuxes between the atmosphere and the ocean are calculated on thesea level grid. In this section, we describe the conversion from the ocean model grid to the sea
surface grid. The standard variables to be converted from the ocean model to the atmosphericsea surface grid are sea surface temperature (SST ), sea ice concentration (AI), sea ice thickness(HI), snow depth over sea ice (HSN ), sea ice internal temperature (T I), and ocean surface
current velocity (U O, V O). In order to clarify which grid we are dealing with in the future,
variables in the ocean model grid will be denoted by superscript OGCM and variables in thesea surface grid by superscript oc. In addition, the position in the ocean grid is denoted by LO
and the position in the sea surface grid by LC. The boundary condition of the ocean in the sea
level grid is deﬁned as follows.
X

IJO(LC)oc

SST (LC) =[SST OGCM (IJO2C(LC, N )) × SOCN (LC, N )]/SOCN G(LC) (4.13)

N =1

X

IJO(LC)oc

AI (LC) =[AI OGCM (IJO2C(LC, N )) × SOCN (LC, N )]/SOCN G(LC)(4.14)N =1

1
HI (LC) =SOCN G(LC) × AI oc (LC)X

IJO(LC)oc

×AI

HSN oc (LC) =OGCM[HI OGCM (IJO2C(LC, N ))N =1

(IJO2C(LC, N )) × SOCN (LC, N )]1
SOCN G(LC) × AI oc (LC)X

IJO(LC)[HSN OGCM (IJO2C(LC, N ))N =1

×AI OGCM (IJO2C(LC, N )) × SOCN (LC, N )]1
T I (LC) =SOCN G(LC) × HI oc (LC) × AI oc (LC)
oc

×HI

OGCM

(IJO2C(LC, N )) × AI
PIJO(LC)U O (LC) = RU O(LC) ×
oc

N =1

PIJO(LC)+RV O(LC) ×

N =1

V O (LC) = −RV O(LC) ×

N =1

N =1

IJO(LC)[T I OGCM (IJO2C(LC, N ))N =1

(IJO2C(LC, N )) × SOCN (LC, N )](4.17)[V OOGCM (IJO2C(LC, N )) × SOCN (LC, N )](4.18)SOCN G(LC)PIJO(LC)+RU O(LC) ×

OGCM

X

(4.16)[U OOGCM (IJO2C(LC, N )) × SOCN (LC, N )]SOCN G(LC)PIJO(LC)oc(4.15)[U OOGCM (IJO2C(LC, N )) × SOCN (LC, N )]SOCN G(LC)[V OOGCM (IJO2C(LC, N )) × SOCN (LC, N )](4.19)SOCN G(LC)126X

IJO(LC)SOCN G(LC) =SOCN (LC, N )(4.20)

N =1

where, IJO(LC): Number of ocean grids corresponding to the sea level grid (LC) in theatmospheric node.
IJO2C(LC, N ): Location of the ocean grid corresponding to the sea surface grid in theatmospheric node.
SOCN (LC, N ): Area of the ocean grid corresponding to the sea surface grid in the atmospheric node.
RU O(LC): Cosine of the rotation angle of the vector.
RV O(LC): Sine of the rotation angle of the vector
SOCN G(LC): Area of ocean occupied by sea surface grid.
The ratio of land surface to the sea level grid is also deﬁned as follow.
F LN Doc = (1 − SOCN G)/SO
where, SO is the area of sea surface grid.
The variables related to sea ice that are converted to the sea surface grid are calculated
as the average of variables categorized by sea ice layer thickness (AIM, HIM, HSM, T IM ) as
follows.

AI OGCM =NICX

AIM OGCM (L)(4.21)HIM OGCM (L) × AIM OGCM (L)/AI OGCM

(4.22)

L=1

HI OGCM =

N
IC
X
L=1

N
IC
X

HSM OGCM (L) × AIM OGCM (L)/AI OGCM

(4.23)T IM OGCM (L) × AIM OGCM (L)/(AI OGCM × HI OGCM )(4.24)HSN OGCM =

L=1

TI

OGCM

=

N
IC
X
L=1

where, N IC is the number of category of sea ice.
4.1.3.2

Conversion of Air-Sea Fluxes Calculated on the Sea Surface Grid to the
Ocean Grid

Fluxes calculated on the sea surface grid are calculated at sea surface and sea ice surface,respectively, and ﬂuxes to the atmosphere are calculated as
F LXO = (1 − AI) × F LU XO + AI × F LU XI(4.25)

These ﬂuxes are time-integrated by the ﬂux coupler in the atmospheric model with weights forsea surface and sea ice extent, and then converted to the ocean grid by the coupled atmosphereocean time step and passed to the ocean model.
127

F LU XOAOGCM (LO) = ROCN (LO)X

IJA(LO)×F LU XOAoc (IJC2O(LO, N )) × SAT M (LO, N )SAT M G(LO)N =1(4.26)F LU XIAOGCM (LO) = ROCN (LO)X

IJA(LO)×

N =1

F LU XIAoc (IJC2O(LO, N )) × SAT M (LO, N )SAT M G(LO)X(4.27)IJA(LO)SAT M G(LO) = ROCN (LO) ×

SAT M (LO, N )(4.28)N =1

ROCN (LO) = SAT M G(LO)/S OGCM (LO)

(4.29)F LU XOAoc = (1 − AI oc ) × F LU XOoc

(4.30)F LU XIAoc = AI oc × F LU XI oc(4.31)

where, IJA(LO)：Number of sea level grids in the atmospheric model corresponding to theocean grid (LO)IJC2O(LO, N )：Location of the sea surface grid of the atmospheric model corresponding
to the ocean grid
SAT M (LO, N )：Area of the sea surface grid of the atmospheric model corresponding to theocean grid
SAT M (LO, N ) = SOCN (LC, L), LC = IJC2O(LO, N ), LO = IJO2C(LC, L)
S OGCM (LO)：Area of the ocean grid
The area of the ocean grid and the sum of the areas of the corresponding ocean grids shouldmatch, although the coordinate systems of the atmospheric model and the ocean model are
diﬀerent (the surface areas of the earth in the atmospheric model and the ocean model do notmatch exactly).
When creating the conversion ﬁle, the grid of the atmospheric model is divided into small
areas, and the area of the corresponding ocean grid is estimated from the sum of these areas,so they do not match exactly. For this reason, the ﬂux balance between the atmosphere andthe ocean is adjusted by multiplying by the ratio (ROCN ). The wind stresses to the ocean are
also calculated as wind stresses over sea level (T XO, T Y O) and over sea ice (T XI, T Y I), but
without multiplying the weights of sea level and sea ice area.
T XOOGCM (LO) =X

IJA(LO)+RU (LO) × ROCN (LO) ×

N =1

X

IJA(LO)+RV (LO) × ROCN (LO) ×

N =1[T XOoc (IJC2O(LO, N )) × SAT M (LO, N )]SAT M G(LO)[T Y Ooc (IJC2O(LO, N )) × SAT M (LO, N )]SAT M G(LO)128(4.32)

T Y OOGCM (LO) =X

IJA(LO)−RV (LO) × ROCN (LO) ×

N =1

X

IJA(LO)+RU (LO) × ROCN (LO) ×

N =1[T XOoc (IJC2O(LO, N )) × SAT M (LO, N )]SAT M G(LO)[T Y Ooc (IJC2O(LO, N )) × SAT M (LO, N )]SAT M G(LO)(4.33)

where,RU (LO)：cosine of the rotation angle of the vector
RV (LO)：sine of the rotation angle of the vector
4.1.3.3

Redistribution of Fluxes in the Ocean Model

The ﬂuxes converted to the ocean grid are updated at each time step of the coupling. Since
the coupling time step is longer than the ocean model time step, the sea level/sea ice area
ratio in the ocean model is updated to a diﬀerent value than the one used to calculate theﬂux. Therefore, in order to obtain an accurate heat and water balance, the ﬂuxes need to be
distributed according to the updated sea surface and sea ice area ratios. The ﬂuxes F LU XOA
and F LU XIA at sea surface and sea ice surface are ocean grid-averaged values. Fluxes to eachsea ice category are currently distributed evenly independent of sea ice thickness.
F LU XIAM (L) = F LU XIA × AIM (L)F LU XOA = F LU XOA + F LU XIA × [1.0 −

LM
AX
X

(4.34)AIM (L)](4.35)

L=1

where, AIM denotes the percentage of area covered by sea ice in the grid (sea ice concentration), L denotes the type of sea ice thickness category, and LM AX denotes the number ofthickness categories. If there is no sea ice surface, all ﬂuxes will be at sea surface. If sea ice
disappears in the middle of the coupling time step, the ﬂux due to sublimation is divided into
the heat ﬂux assuming a sea ice surface and the freshwater ﬂux (sea ice loss). The heat ﬂux is
directly reﬂected in the temperature change of the ﬁrst layer of ocean. On the other hand, thefreshwater ﬂux due to sublimation is converted into heat ﬂux and freshwater ﬂux and given to
the ﬁrst layer of the ocean, assuming that sea ice is generated by the sublimation. As for thewind stress, it is not weighted by sea level and sea ice area before the grid transformation, so it
is driven by the respective area weights in each sea ice thickness category in the ocean model.
For this reason, momentum is not conserved.
4.1.3.4

Water Runoﬀ from Rivers to the Ocean

At the end of the river model, we calculate the water ﬂowing from the estuary of river tothe ocean. Water arriving at the estuary of the river grid is ﬁrst converted to the atmosphericsea surface grid and time integrated in a ﬂux coupler. After that, it is converted to the ocean
grid via an exchanger and passed to the ocean model in the same way as the atmosphericprecipitation data. At this point, the temperature of the river water is treated as the same as
the sea surface temperature, as is the case with precipitation. Therefore, strictly speaking, heatis not conserved. Ice runoﬀ is handled in the same way as snowfall.
1294.1.3.5

Number of Divisions in the Sea Surface Grid and Resolution of the Ocean
Model

The sea surface grid is created by dividing the latitude and longitude of the atmospheric grid,but if the number of divisions is not suﬃcient and the ocean model grid has a higher resolution
than the atmospheric sea surface grid, the structure of the atmospheric grid size may remain
when the ﬂux is converted to the ocean grid through the exchanger. . In addition, data suchas precipitation from the atmosphere is not interpolated when converting from the atmosphericgrid to the ocean grid, so the atmospheric grid structure remains in the ocean grid for theseﬂuxes. When linear interpolation is used instead of cubic spline interpolation when convertingto the sea surface grid, the atmospheric grid structure may remain for diﬀerential quantities
such as wind stress curl.

130

4.2

Deﬁnition of Land-Sea Distribution

The land-sea distribution in MIROC is prioritized by the land-sea distribution deﬁned by theocean model. While one grid in the ocean model is deﬁned by land or sea only, the land andocean grids in the atmospheric model are determined in proportion to the land and sea to beconsistent with the ocean model’s land-sea distribution.
SA : area of the atmospheric grid, SLij : area of the land grid, SOij : area of the sea surface
land , F LN D oc : percentage of land surface is occupied by each grid.
grid, F LN Datm , F LN Dij
ij
Then, following equation is satisﬁed.

SA ∗ F LN D

atm

=

jldiv
X ildiv
X

land(SLij ∗ F LN Dij)=

j=1 i=1

jodiv
X iodiv
X

oc(SOij ∗ F LN Dij)(4.36)

j=1 i=1

where, (ildiv,jldiv) is the number of east-west and north-south divisions of the land surface
grid, and (iodiv,jodiv) is the number of east-west and north-south divisions of the sea surface
grid. In the land surface grid, if even a small amount of land is deﬁned to exist, boundary values
such as land cover are required.

131

5

References
1. Abdul-Razzak, H., and Ghan, S. J. (2000). A parameterization of aerosol activation: 2.
Multiple aerosol types. Journal of Geophysical Research: Atmospheres.
https://doi.org/10.1029/1999jd901161
2. Arakawa, A., , and W. H. Schubert, (1974). Interaction of a cumulus cloud ensemble withthe large-scale environment, Part I. Journal of Atmospheric Science, 31 , 674–701.
3. Arakawa, A., , and K.-M. Xu, (1990). The macroscopic behavior of simulated cumulus
convection and semiprognostic tests of the Arakawa–Schubert cumulus parameterization.
Proc. Indo-U.S. Seminar on Parameterization of Subgrid-Scale Processes in Dynamical
Models of Medium-Range Prediction and Global Climate, Pune, India, Indian Institute ofTropical Meteorology, 3–18.
4. Arakawa, A., and K.-M. Xu, (1992). The macroscopic behavior of simulated cumulus
convection and semiprognostic tests of the Arakawa-Schubert cumulus parameterization.
Physical Process in Atmospheric Models., 3–18.
5. Arakawa, A., and Konor, C. S. (1996). Vertical Diﬀerencing of the Primitive Equations
Based on the Charney–Phillips Grid in Hybrid andsigma–p Vertical Coordinates. Monthly
Weather Review, 124(3), 511–528.
6. Asselin, R. (1972). Frequency Filter for Time Integrations. Monthly Weather Review,100(6), 487–490.
7. Berry, E. X. (1968). Modiﬁcation of the warm rain process. In Proc. First Conf. on
Weather Modiﬁcation, AMS, 1968, 81-85.
8. Berger, A. (1978). Long-term variations of caloric insolation resulting from the Earth’s
orbital elements. Quaternary research, 9(2), 139-167.
9. Bourke, W. (1988). Spectral Methods in Global Climate and Weather Prediction Models.
In M. E. Schlesinger (Ed.), Physically-Based Modelling and Simulation of Climate andClimatic Change: Part 1 (pp. 169–220). Dordrecht: Springer Netherlands.

10. Bretherton, C. S., McCaa, J. R., and Grenier, H. (2004). A new parameterization for shallow cumulus convection and its application to marine subtropical cloud-topped boundary
layers. Part I: Description and 1D results. Monthly Weather Review, 132(4), 864–882.
11. Bushell, A. C., Wilson, D. R., and Gregory, D. (2003). A description of cloud production
by non-uniformly distributed processes. Quarterly Journal of the Royal MeteorologicalSociety, 129(590), 1435–1455.
12. Businger, J. A., Wyngaard, J. C., Izumi, Y., and Bradley, E. F. (1971). Flux-ProﬁleRelationships in the Atmospheric Surface Layer. Journal of the Atmospheric Sciences,28(2), 181–189.
13. Cesana, G., Waliser, D. E., Jiang, X., and Li, J. ‐L F. (2015). Multimodel evaluation ofcloud phase transition using satellite and reanalysis data. Journal of Geophysical Research,
120(15), 7871–7892.
14. Chikira, M., and Sugiyama, M. (2010). A Cumulus Parameterization with State-Dependent
Entrainment Rate. Part I: Description and Sensitivity to Temperature and Humidity Proﬁles. Journal of the Atmospheric Sciences, 67(7), 2171–2193.
15. Colella, P., and Woodward, P. R. (1984). The Piecewise Parabolic Method (PPM) for
gas-dynamical simulations. Journal of Computational Physics.

132https://doi.org/10.1016/0021-9991(84)90143-8
16. Diehl, K., Simmel, M., and Wurzler, S. (2006). Numerical sensitivity studies on theimpact of aerosol properties and drop freezing modes on the glaciation, microphysics, anddynamics of clouds. Journal of Geophysical Research, 111(D7).
https://doi.org/10.1029/2005jd005884
17. Gregory, D. (2001). Estimation of entrainment rate in simple models of convective clouds.
Quarterly Journal of the Royal Meteorological Society, 127(571), 53–72.
18. Gregory, D., Kershaw, R., and Inness, P. M. (1997). Parametrization of momentum transport by convection. II: Tests in single-column and general circulation models. Quarterly
Journal of the Royal Meteorological Society, 123(541), 1153–1183.
19. Haltiner, G. J., and Williams, R. T. (1980). Numerical prediction and dynamic meteorology (No. 551.5 HAL ). sidalc.net. Retrieved fromhttp://www.sidalc.net/cgi-bin/wxis.exe/?IsisScript=FCL.xisandmethod=
postandformato=2andcantidad=1andexpresion=mfn=004222
20. Hasumi, H. (2015). CCSR ocean component model (COCO) (Version 4.0) (p. 68). Atmosphere and Ocean Research Institute, The University of Tokyo. Retrieved fromhttps://ccsr.aori.u-tokyo.ac.jp/ hasumi/COCO/coco4.pdf
21. Helfand, H. M., and Labraga, J. C. (1988). Design of a Nonsingular Level 2.5 Second-Order
Closure Model for the Prediction of Atmospheric Turbulence. Journal of the Atmospheric
Sciences, 45(2), 113–132.
22. Hines, C. O. (1997) Doppler-spread parameterization of gravity wave momentum deposition in the middle atmosphere, Part 2: Broad and quasi monochromatic spectra, andimplementation, Journal of Atmosphere Sol.-Terr. Physics, 59, 387–400.
23. Holtslag, A. A. M., and Boville, B. A. (1993). Local Versus Nonlocal Boundary-Layer
Diﬀusion in a Global Climate Model. Journal of Climate, 6(10), 1825–1842.
24. K-1 model developers (2004). K-1 coupled GCM (MIROC) description, K-1 Tech. Rep.,
1, edited by: Hasumi, H. and Emori, S., Center for Climate System Research, the Univ.
of Tokyo, Tokyo, 34 pp.
25. Kain, J. S., and Michael Fritsch, J. (1990). A One-Dimensional Entraining/Detraining
Plume Model and Its Application in Convective Parameterization. Journal of the Atmospheric Sciences, 47(23), 2784–2802.
26. Lew, J. K., , D. C. Montague, , and H. R. Pruppacher, (1986) A wind tunnel investigationon the riming of snowﬂakes. Part I: Porous disks and large stellars. Journal of Atmospheric
Science, 43, 2392–2409.
27. Lin, S.-J., and Rood, R. B. (1996). Multidimensional Flux-Form Semi-Lagrangian Transport Schemes. Monthly Weather Review, 124(9), 2046–2070.
28. Lohmann, U. (2004). Can anthropogenic aerosols decrease the snowfall rate?. Journal of
the atmospheric sciences, 61(20), 2457-2468.
29. Lohmann, U., and Diehl, K. (2006). Sensitivity studies of the importance of dust ice
nuclei for the indirect aerosol eﬀect on stratiform mixed-phase clouds. Journal of the
Atmospheric Sciences, 63(3), 968–982.
30. Louis, J., Tiedtke, M., and Geleyn, J. (1982). A short history of the PBL parameterization
at ECMWF, paper presented at Workshop on Planetary Boundary Layer Parameterization, Eur. Cent. for Medium Range Weather Forecasts. Reading, UK.
133

31. Louis, J.-F. (1979). A parametric model of vertical eddy ﬂuxes in the atmosphere. BoundaryLayer Meteorology, 17(2), 187–202.
32. McFarlane, N. (1987). The Eﬀect of Orographically Excited Gravity Wave Drag on the
General Circulation of the Lower Stratosphere and Troposphere. Journal of the Atmospheric Sciences, 44(14), 1775-1800.
33. Mellor, G. L. (1973). Analytic Prediction of the Properties of Stratiﬁed Planetary Surface
Layers. Journal of the Atmospheric Sciences, 30(6), 1061–1069.
34. Mellor, G. L., and Yamada, T. (1974). A Hierarchy of Turbulence Closure Models forPlanetary Boundary Layers. Journal of the Atmospheric Sciences, 31(7), 1791–1806.
35. Mellor, G. L., and Yamada, T. (1982). Development of a turbulence closure model for
geophysical ﬂuid problems. Reviews of Geophysics, 20(4), 851-875.
36. Mesinger, F., and Arakawa, A. (1976). Numerical methods used in atmospheric models(GARP Publications Series) (p. 64). CAU. Retrieved fromhttp://eprints.uni-kiel.de/40278/
37. Miller, M.J., A.C.M. Beljaars and T.N. Palmer, (1992). The sensitivity of the ECMWF
model to the parameterization of evaporation from the tropical oceans. Journal of Climate,5, 418-434.
38. Miura, H. (2002). Vertical diﬀerencing of the primitive equations in a σ − p hybrid
coordinate (For Spectral AGCM)) (Version DRAFT). Center for Climate System Research,University of Tokyo.
39. Nakajima, T., Tsukamoto, M., Tsushima, Y., Numaguti, A., and Kimura, T. (2000).
Modeling of the radiative process in an atmospheric general circulation model. Applied
Optics, 39(27), 4869–4878.
40. Nakanishi, M. (2001). IMPROVEMENT OF THE MELLOR–YAMADA TURBULENCE
CLOSURE MODEL BASED ON LARGE-EDDY SIMULATION DATA. Boundary-Layer
Meteorology, 99, 349–378.
41. Nakanishi, M., , and H. Niino, (2004). An improved Mellor–Yamada level-3 model withcondensation physics: Its design and veriﬁcation. Boundary-Layer Meteor, 112 , 1–31.
42. Nitta, T., Yoshimura, K., Takata, K., O ishi, R., Sueyoshi, T., Kanae, S., et al. (2014).
Representing Variability in Subgrid Snow Cover and Snow Depth in a Global Land Model:
Oﬄine Validation. Journal of Climate, 27(9), 3318–3330.
43. Ogura, T. (2015). Implementation of shallow cumulus Parameterization. National Institute of Environmental Studies.
44. Ogura, T., Shiogama, H., Watanabe, M., Yoshimori, M., Yokohata, T., Annan, J. D., Hargreaves, J. C., Ushigami, N., Hirota, K., Someya, Y., Kamae, Y., Tatebe, H., and Kimoto,M. (2017). Eﬀectiveness and limitations of parameter tuning in reducing biases of topof-atmosphere radiation and clouds in MIROC version 5, Geoscience Model Development,
10, 4647–4664
45. Pan, D.-M. (1995, January 1). Development and Application of a Prognostic Cumulus
Parametrization. ui.adsabs.harvard.edu. Retrieved fromhttps://ui.adsabs.harvard.edu/abs/1995PhDT........74P
46. Pan, D-M., , and D. A. Randall, 1998: A cumulus parameterization with a prognostic
closure. Quarterly Journal of the Royal Meteorological Society, 124 , 949–981.

13447. Park, S., and Bretherton, C. S. (2009). The University of Washington Shallow Convection and Moist Turbulence Schemes and Their Impact on Climate Simulations with theCommunity Atmosphere Model. Journal of Climate, 22(12), 3449–3469.
48. Pruppacher, H. R., , and J. D. Klett, (1997). Microphysics of Clouds and Precipitation.
Kluwer Academic, 954 pp.
49. Randall, D. A. (1980). Conditional Instability of the First Kind Upside-Down. Journal of
the Atmospheric Sciences, 37(1), 125–130.
50. Randall, D. A., and Pan, D.-M. (1993). Implementation of the Arakawa-Schubert Cumulus
Parameterization with a Prognostic Closure. In K. A. Emanuel and D. J. Raymond (Eds.),
The Representation of Cumulus Convection in Numerical Models (pp. 137–144). Boston,
MA: American Meteorological Society.
51. Randall, D. A., , D-M. Pan, , P. Ding, , and D. G. Cripe, (1997). Quasi-equilibrium. The
Physics and Parameterization of Moist Atmospheric Convection, R. K. Smith, Ed., Kluwer
Academic, 359–386.
52. Sekiguchi, M., and Nakajima, T. (2008). A k-distribution-based radiation code and its
computational optimization for an atmospheric general circulation model. Journal of
Quantitative Spectroscopy and Radiative Transfer, 109(17), 2779–2793. vSmith, R. K.
(2013). The Physics and Parameterization of Moist Atmospheric Convection. Springer
Science and Business Media.
53. Sommeria, G., and Deardorﬀ, J. W. (1977). Subgrid-Scale Condensation in Models ofNonprecipitating Clouds. Journal of the Atmospheric Sciences, 34(2), 344–355.
54. Stevens, B. (2005). ATMOSPHERIC MOIST CONVECTION. Annual Review of Earthand Planetary Sciences, 33(1), 605–643.
55. Suzuki, T., Saito, F., Nishimura, T., and Ogochi, K. (2009). Coupling procedures of heatand freshwater ﬂuxes in the MIROC (Model for Interdisciplinary Research on Climate)version 4. JAMSTEC report of research and development
56. Takata, K., Emori, S., and Watanabe, T. (2003). Development of the minimal advanced
treatments of surface interaction and runoﬀ. Global and Planetary Change, 38(1), 209–222.
57. Takemura, T. (2005). Simulation of climate response to aerosol direct and indirect eﬀects
with aerosol transport-radiation model. Journal of Geophysical Research, 110(D2).
https://doi.org/10.1029/2004jd005029
58. Takemura, T., Okamoto, H., Maruyama, Y., Numaguti, A., Higurashi, A., and Nakajima,
T. (2000). Global three-dimensional simulation of aerosol optical thickness distribution ofvarious origins. Journal of Geophysical Research, 105(D14), 17853–17873.
59. Takemura, T., Nakajima, T., Dubovik, O., Holben, B. N., and Kinne, S. (2002). Singlescattering albedo and radiative forcing of various aerosol species with a global threedimensional model. Journal of Climate, 15(4), 333–352.
60. Takemura, T., Egashira, M., Matsuzawa, K., Ichijo, H., O ishi, R., and Abe-Ouchi, A.
(2009). A simulation of the global distribution and radiative forcing of soil dust aerosols
at the Last Glacial Maximum. Atmospheric Chemistry and Physics, 9(9), 3061–3073.
61. Tatebe, H., Ogura, T., Nitta, T., Komuro, Y., Ogochi, K., Takemura, T., et al. (2019).
Description and basic evaluation of simulated mean state, internal variability, and climate
sensitivity in MIROC6. Geoscientiﬁc Model Development, 12(7), 2727–2765.
62. Tompkins, A. M. (2002). A prognostic parameterization for the subgrid-scale variability of135

water vapor and clouds in large-scale models and its use to diagnose cloud cover. Journalof the Atmospheric Sciences, 59(12), 1917–1942.
63. Watanabe, M., Emori, S., Satoh, M., and Miura, H. (2009). A PDF-based hybrid prognostic cloud scheme for general circulation models. Climate Dynamics, 33(6), 795–816.
64. Watanabe, S. (2008). Constraints on a Non-orographic Gravity Wave Drag Parameterization Using a Gravity Wave Resolving General Circulation Model, SOLA, 4, 61-64.
Watanabe, S., Hajima, T., Sudo, K., Nagashima, T., Takemura, T., Okajima, H., Nozawa,T., Kawase, H., Abe, M., Yokohata, T., Ise, T., Sato, H., Kato, E., Takata, K., Emori,S., and Kawamiya, M. (2011) MIROC-ESM 2010: model description and basic results of
CMIP5-20c3m experiments, Geoscience Model Development, 4, 845-872.
65. Williams, P. D. (2009). A Proposed Modiﬁcation to the Robert–Asselin Time Filter.
Monthly Weather Review, 137(8), 2538–2546.
66. Wilson, D. R., and Ballard, S. P. (1999). A microphysically based precipitation scheme for
the UK meteorological oﬃce uniﬁed model. Quarterly Journal of the Royal MeteorologicalSociety, 125(557), 1607-1636.
67. Wood, R., and Bretherton, C. S. (2006). On the relationship between stratiform low cloud
cover and lower-tropospheric stability. Journal of Climate, 19(24), 6425-6432.
68. Xu, K.-M. (1991). The coupling of cumulus convection with large-scale processes (Ph.D.).
University of California. Retrieved fromhttps://elibrary.ru/item.asp?id=5851863
69. Xu, K.-M., (1993). Cumulus ensemble simulation. The Representation of Cumulus Convection in Numerical Models, Meteor. Monogr., American Meteorological Society,46, 221235.

136