NICAM lecture
Masaki Satoh2012/02/23

contents• Dynamical core– governing equations– Time integration• Horizontal grid structure• Vertical grid structure

DYNAMICAL CORE

Non‐hydrostatic equations
 Flux form equations
 mass momentum
 Total energy

 No approximationq     qv   St Fully compreible nonhydrostatic equations
 Deep atmosphere Time integration
 Split-explicit method
 Small-step: horizontally explicit, vertically implicit
 Large-step: 2nd, 3rd-order Runge-Kutta method

Governing equations

  W


3
R  h  Vh  1/ 2  G  Vh   0t  GVh  h P 
G 3 P  ADVh  FCoriolist  P W
 1/ 2   Rg  ADVz  Fz ,Coriolist  G 






 
 W

3
h  k    1 / 2  G  Vh   0
Etotal  h  h  k   Vh  

t
 

G
Prognostic variables• Density• Horizontal momentum• Radial momentum• (Internal energy)• Total energy

RG 
1/ 2

Vh  G1/ 2  vh
W  G1 / 2  w
E  G1/ 2  eininternal energy + kinetic energy + potential energy
Etotal   G 1 / 2 ein  k   

Metrics

 z 
G1/ 2      x , y
G 3  h zH ( z  zs )H  zsTemporal scheme

1.
2.
3.
4.

When values at t=A are known, the slow modetendency S(A) is calculated.
Integrate from A to B using S(A) with fast modetime step.
The slow mode tendency at t=B is calclated:
S(B).
Then, integrate from A to C using S(B) with fastmode time step.

Slow mode: 2nd order Runge‐Kutta or
3rd order Runge‐Kutta
Fast mode: Euler method&TIMEPARAMINTEG_TYPE = 'RK2',
LSTEP_MAX = 28800,
SSTEP_MAX = 4,
DTL = 30.0D0,
SPLIT=.TRUE.,
start_year = 2009,
start_month = 7,
start_day = 1,/

Split explicit scheme: Small Step integration 1
3steps for the time splitting scheme
1. Horizontal Explicit Step
 Horizontal momentum 水平運動量の更新

2. Vertical Implicit Step
 Vertical momentum and density 鉛直運動量と密度の更新

3. Energy Correction Step
 Total energy 内部エネルギーの更新

Horizontal Explicit Step
– Horizontal momentum
t  n
[ t , or t  t / 2 ]



V


h


3
t  ( n 1) t  n
G P    t Vh
 Vh    h P slow mode [ t , or t  t / 2 ]Where Vh  t slow modeis the Slow Mode Tendency

Split explicit scheme: Small Step integration 2
Vertical Implicit Step
– Using horizontal momentum obtained by the Horizontal Explicit Step,
Eqs. Of R,W,E can be written as  W t  ( n 1)  
R t  ( n 1)   R t  n  GR
1/ 2
  G  P t  ( n 1)  
W t  ( n 1)   W t  nt  ( n 1) 




R
g

  G1/ 2 

 Gz

  W t  ( n 1)   2 t  n  Rd t  ( n 1)  ~P t  ( n 1)   P t  n
cs

W
g 


1/ 2

  G

 CV

Rd
GE
CV

(6)
(7)
(8)

– The Helmholtz equation for w can be made by combining Eqs. (6)‐(9)：  W t  ( n 1)  
  1   2 2 t  n W t  ( n 1)      2 Rd ~ W t  ( n 1)  2   g      csg
G1/ 2    CV
G1/ 2   G1/ 2   G1/ 2   W t  n   Gz  Rd
  1   t  n
t  n


P



G
  GR g
 1/ 2E     R  G  
CV






(9)

– Eq. (9) is soleved for W– Eq. (6) gives R– Instead of Eq. (8), energy E is solve by the total energy equation, next.

Split explicit scheme: Small Step integration 3
Energy Correction Step
(Total energy) = (Internal energy) + (kinetic energy) + (potential energy)– Using the conservation law of the total energy equation

 

 W
3


Etotal  h  h  k   Vh  
h
k
G
V
h
 0 (10)








1/ 2
t
 
G

Etotal   G 1 / 2 ein  k   
where– Eq.(10) is discretized as follows

 
 W
 
t  ( n 1) t  n
3
EtotalhkGVh  Etotal h  h  k   Vh  








1/ 2

G







t  ( n 1) This gurantees theconservation of the total energy– Since Kinetic energy and potential energy are alreadyknown, the internal energy is obtained, and thus
temperature and pressure can be obtained.

t  ( n 1) 
E t  ( n 1)   Etotal  t  ( n 1)  G 1 / 2 k t  ( n 1)   Time integration of the Large Step
Two steps of the Large time step integration
Coriolis terms
1. Advection terms– Calculation of Advection of momentum– 3dim velocity V is seprated into
horizontal velocity Vh and verticalvelocity W– Vh is expressed by the Cartesiancomponents fixed to the earth center(V1,V2,V3)– Divergence of each components  V1 ,   V2 ,   V3 
– Summing up the horizontal and verticalcomponents of Advection ADVh , ADVz

Expression of momentum• Use 4 variables of momentum (u , v, w)

Horizontal velocity(u , v)

Vertical velocityw

Three components of the
Cartesian coordinate are used torepresent the horizontal velocity.

(v x , v y , v z )HORIZONTAL GRID

Horizontal grid division: glevel

Number of Grid points

N grid  10 4

glevelNumber of Grid points
Definition of horizontal resolution

N grid  10 4
g‐levelglevel1/ 2

 4R x  
N 
 grid 
2E

Grid Number

Horizontal Resolution [km]5

10,240223.2

640,960111.6

7163,84055.80

8655,36027.90

92,621,440

13.95

10

10,485,760

6.974

11

41,943,040

3.487

12

167,772,160

1.744

13

671,088,640

0.872

14

2,684,354,560

0.440

Discretization of sphere

Control volume
Hexagon or pentagon

Definition ofdivergence operator

Problem of original icosahedral grid
Error of divergence operator(rigid rotation field)•

Large error along the arcs of theoriginal icosahedron• Fractal distribution• Sources of grid noise

Distribution of area of control volume

We expect that smoothing of CV
areas will reduce errors.

Modification of icosahedral grid
1. Spring dynamics smoothing
6
dw 0
M  k (d i  d )e i  w 0
dt
i 1
2a
dr0
 w0
d 
dt
10  2l 1

mkgrid2. Change of definition point toareal center
Second order accuracy

mkgcgridOriginal icosahedral grid Modified icosahedral grid

Tomita, Satoh, Goto (2002, J. Comp. Phys.)Horizontal uniformity vs 
格子間隔の最大最小比 lmax / lmin

In mkgrid.cnf
spr_beta=1.15D0Standard grid（ＳＴＤ-grid）
lmax/lmin  1.34.

Spring grid（SPR-grid）
lmax/lmin ：  にも glevelにも依存If  =1.2, lmax/lmin  1.24.

Area distributions of control volume
Original icosahedral gridMdified grid

Tomita et al. (2001, J. Comp. Phys.)

PARALLEL PROCESSING

Parallel strategy(0) rlevel-0(1) rlevel-1(2) rlevel-2(3) rlevel-3

R‐level

ADM_lall  10  4

ADM_gall 

rlevel

N grid
ADM_lall.

 4( glevel rlevel ).

Concept of Horizontal Grid
• What is “rlevel” ?
– Region Division Level– Designed for MPI processes

rlevel 0
10 regions / globe

rlevel 1
40 regions / globe

rlevel # of regions
0
10
1
40
2
160
3
640
4
2560
5
10240

If using 20 processes,
128 regions / process– File name: var‐name.rgn????? e.g. sst.rgn00000 ~ sst.rgn00639– Note: As rlevel increases, total memory consumption is increased, especiallywhen glevel‐rlevel is small.
• due to halo (overlapped) region.

Exmple: Rlevel‐1

Indexing of the horizontal grid
 Region structure: diamonds
Structured grid
 Fortran 2D array
Efficient for the vector machine 2d arrays are onedimensionalized to gain longer
vector length.
Vector ratio : more than ９９％

Typical Arrays
V(ij,k,l)Q(ij,k,I,nq)

Grid indexing

ADM_gall_1d  2 ( glevel  rlevel )  2,



ADM_gall  ADM_gall_1d   2
2( glevel rlevel )

ADM_gall_in  ADM_gall_1d  1  2
2

ADM_gmin  1  1,
ADM_gmax  1  ( 2 ( glevel rlevel ) ).

( glevel rlevel )



2

2 .



2

1 .

Prognostic variables! Prognostic variables!
!‐‐‐ rho X ( G^{1/2} X gamma2 )Real(8) :: rhog(ADM_gall_in,ADM_kall,ADM_lall)!‐‐‐ rho X ( G^{1/2} X gamma2 ) X vx
Real(8) :: rhogvx(ADM_gall_in,ADM_kall,ADM_lall)!‐‐‐ rho X ( G^{1/2} X gamma2 ) X vy
Real(8) :: rhogvy(ADM_gall_in,ADM_kall,ADM_lall)!‐‐‐ rho X ( G^{1/2} X gamma2 ) X vz
Real(8) :: rhogvz(ADM_gall_in,ADM_kall,ADM_lall)!‐‐‐ rho X ( G^{1/2} X gamma2 ) X w
Real(8) :: rhogw(ADM_gall_in,ADM_kall,ADM_lall)!‐‐‐ rho X ( G^{1/2} X gamma2 ) X ein
Real(8) :: rhoge(ADM_gall_in,ADM_kall,ADM_lall)!‐‐‐ rho X ( G^{1/2} X gamma2 ) X eintot
Real(8) :: rhogetot(ADM_gall_in,ADM_kall,ADM_lall)!‐‐‐ rho X ( G^{1/2} X gamma2 ) X q
Real(8) :: rhogq(ADM_gall_in,ADM_kall,ADM_lall,TRC_VMAX)

Vertical grid
Terrain following coordinate in general form   z
Only the following simplest form is assumed now.

zT  z  z s zT  z s  zT   zT ;   z s   0
zT  z sz  zs zT

GRD_gzh ( ADM_kall)  GRD_htop
GRD_gzh ( ADM_kmin )  0

Vertical layer depth:
GRD_dgz (k )  GRD_gzh (k  1)  GRD_gzh (k ),GRD_dgzh (k )  GRD_gz (k )  GRD_gz (k  1),
GRD_dgz (ADM_kall)  GRD_dgz(ADM_kall  1),
GRD_dgzh (1)  GRD_dgzh (ADM_kmin ).

Conversion from real height to model height:
GRD_vz(ij,k,l_region, GRD_Z)  GRD_zs(ij,l_region)
 GRD_htop,
GRD_htop  GRD_zs(ij,l_region)
GRD_vz(ij , k , l _ region, GRD_ZH)  GRD_zs(ij,l_region)
GRD_gzh (k )  GRD_htop,
GRD_htop  GRD_zs(ij,l_region)

GRD_gz(k ) Conversion from model height to real height:
 GRD_htop  GRD_zs
GRD_vz ( k , GRD_Z )  GRD_zs  
GRD_htop  GRD_gz ( k ) GRD_htop - GRD_zs
GRD_vz( k, GRD_ZH)  GRD_zs  
GRD_htop  GRD_gzh ( k )Default vertical grid levels

NICAM database: horizontal grid data
Common data base is prepared:

NICAM_DATABASE/hgrid/glxx/rlyy/grid.rgnZZZZZ
do l_region=1, ADM_lall
write (wc(1:5),’(i5.5)’) l_region‐1 ! region id ZZZZZ

xx: glevel‐xx
yy: rlevel‐yy
ZZZZZ: region id

open (fid, file=”grid.rgn”//trim(wc(1:5)), form=”unformatted”, access=”sequential”)
read (fid) gall_1d! Grid number of of one side of region inc. halo (integer)gall = gall_1d*gall_1d! Total number of region inc. halo (integer)

read (fid) GRD_xt (1:gall,1,l_region,GRD_XDIR) ! X‐cord. of grid point (real8)
read (fid) GRD_xt (1:gall,1,l_region,GRD_YDIR) ! Y‐cord. of grid point(real8)
read (fid) GRD_xt (1:gall,1,l_region,GRD_ZDIR) ! Z‐cord. of grid point (real8)
close (fid)end doVertical grid data
NICAM_DATABASE/vgrid/vgrid40.dat
open(fid, file=”vgrid40.dat”, form=”unformatted”,
access=”sequential”)
read(fid) vlayer! Total number of vertical layer (integer)kall = vlayer + 2! Total number of vertical index inc. top andbottom (integer)
allocate ( gz(1:kall) )allocate ( gzh(1:kall) )read (fid) gz(:)! Real height [m] of integer levels (real8)
read (fid) gzh(:)! Real heght [m] of half‐integer levels (real8)
close (fid)Information of NICAM• http://nicam.jp/• NICAM CVS server: nicam.jp• NICAM Mail list
– admin@nicam.jp: NICAM administrators• M. Satoh, H. Tomita+H. Yashiro+H. Yamaura (AICS), H.
Taniguchi(IPRC)– develop@nicam.jp: core developers– contrib@nicam.jp: access to the latest code
– users@nicam.jp: access to the tutorial package• All the lecturers are registered to this ML.

Research collaborations• Can use the NICAM code
• Can use the NICAM data, especially producedwith the Earth Simulator, Athena Project, andK‐computer• Can run NICAM on ES or K

TERMS AND CONDITIONS– Do not distribute the NICAM codes to any third party.
– Please report any bugs you have found.
– Please give us a feedback on the usability of NICAM so as toimprove it.
– Please let us know the achievements you made with NICAM,
so that we can make a list on the NICAM website(http://www.nicam.jp/hiki/?NICAM+Papers)– When you present or publish the results obtained with
NICAM, please cite the publications of NICAM. The relevantpublications for the NICAM researches can be found in thecodes or on the NICAM website.
– If a NICAM team member is heavily involved with yourproject, please acknowledge appropriate members or
consider adding him or her as a co‐author.

NICAM tutorial package
Global NICAM package
http://nicam.jp/msn3631/global_nicam/
id:passwd:Regional NICAM package
http://nicam.jp/msn3609/reg_nicam/
id:passwd:Common data set
http://nicam.jp/msn3609/common_data/id:passwd:Global model&NUMFILTERPARAM
dep_hgrid= .false.,
lap_order_divdamp = 2,
hdiff_type= 'NONLINEAR1',cfact= 2.0D0,
Kh_coef_maxlim = 10.0D+16,
Kh_coef_minlim = 1.0D+16,&TIMEPARAM
ZD_hdiff_nl= 20000D0,
INTEG_TYPE = 'RK3',
divdamp_type = 'DIRECT',
LSTEP_MAX = 216,
alpha_d= 1.0D+16,
SSTEP_MAX = 6,alpha_dv= 0.00D0,
DTL= 1200.0D0,
gamma_v= 0.0D0,
SPLIT= .TRUE.,
gamma_h_lap1 = 10.0D+6,
start_year = 2004,
ZD_hdiff_lap1 = 20000.0D0,
start_month = 6,
hdiff_fact_rho = 0.01D0,
start_day = 1,
hdiff_fact_q = 0.0d0,/
ZD= 25000.0d0,##### <‐‐‐ default : RK2 with SSTEP_MAX = 4;alpha_r= 0.0D0,##### <‐‐‐ But, if unstable, use RK3 with
DEEP_EFFECT= .false.,
SSTEP_MAX=6.
/• Glevel 5: dx=220km

Stretch model (original)• Glevel 5: dx_min=22km&TIMEPARAMINTEG_TYPE = 'RK2',
LSTEP_MAX = 10080,
SSTEP_MAX = 8,
DTL = 60.0D0,
SPLIT=.TRUE.,
start_year = 2009,
start_month = 07,
start_day = 01,
start_hour = 00,/! LSTEP_MAX*DTL / 3600 / 24 = 7 (days)!#### <‐‐‐ default : RK2 with SSTEP_MAX = 4;!#### <‐‐‐ But, if unstable, use RK3 with
SSTEP_MAX=6.

&NUMFILTERPARAM
dep_hgrid=.true.
lap_order_divdamp = 2,
hdiff_type = 'NONDIM_COEF',
gamma_h = 0.0217014,
gamma_h_zdep = 0.0325521D0,
ZD_hdiff = 20000.0D0,
divdamp_type = 'NONDIM_COEF',
alpha_d = 0.0217014,
alpha_dv = 0.00D0,
gamma_v = 0.0217014,
gamma_h_lap1 = 0.0D0,
ZD_hdiff_lap1 = 20000.0D0,
hdiff_fact_rho=0.01D0,
hdiff_fact_q=1.0d0,
ZD = 30000.0d0,
alpha_r = 0.0D0,
DEEP_EFFECT = .false.,/

Stretch model (present)• Glevel 5: dx_min=22km&TIMEPARAMINTEG_TYPE = 'RK3',SSTEP_MAX = 6,
LSTEP_MAX = 1440, ! 1day:1440,
2days:2880, 3days: 4230, 7days:10080,@60sec.
DTL= 15.0D0,
SPLIT= .true.,start_year = 2009,
start_month = 07,
start_day = 01,
start_hour = 00,/!!! GL05/trans_alpha=100/DTL=60sec&NUMFILTERPARAM
dep_hgrid= .true.,
lap_order_divdamp = 2,
hdiff_type= 'NONDIM_COEF',
gamma_h= 0.0217014,
gamma_h_zdep = 0.0325521D0,
ZD_hdiff= 20000.0D0,
divdamp_type = 'NONDIM_COEF',
alpha_d= 0.0217014,
alpha_dv= 0.00D0,
gamma_v= 0.0D0,
gamma_h_lap1 = 0.0D0,
ZD_hdiff_lap1 = 20000.0D0,
hdiff_fact_rho = 0.01D0,
hdiff_fact_q = 1.0d0,
ZD= 30000.0d0,
alpha_r= 0.0D0,
DEEP_EFFECT= .false.,/Numerical diffusion
Horizontal numerical diffusion coefficient KH [m4/s] is forth‐order (i.e. KH ∆2 X)
Nondimensional diffusion coefficient gamma is defined as
KH = gamma dx^4/dt
where dx is the grid interval, and dt is the time step.
The corresponding diffusion time istau = dx4/KH = dt/gamma.
The empirical rule for the diffusion coefficient:
As dx => dx/2, tau => tau/2
This means the diffusion time is proportional to the scale of the meteorologicalphenomena.
Under this condition, if dx => dx/2, dt => dt/2, then gamma should not bechanged.
In other cases such as dx => dx/2, dt => dt‘, thentau = dt/gamma => tau/2 = dt'/gamma'
gamma => gamma' = gamma * ( 2 dt'/dt )Transformation method
Transform function: Schmidt transformation  (1  sin  ) /(1  sin  )  1f ( )  sin 1 
  (1  sin  ) /(1  sin  )  1Parameter 

SchmidtRation dx_max/dx_min– =100• Minimum grid interval: 1/10
• Maximum grid interval: 10• Ration dx_max/dx_min = 100– =16• Minimum grid interval: 1/4• Maximum grid interval: 4• Ration dx_max/dx_min = 16

Transgrid&TRANSGRIDPARAM
input_base = './_grid',
output_base = './_s10grid',
type_of_mapping= 'NON_DST',
trans_alpha = 100.0D0,/

Examples of streched grid
Top view

=100

=1000Side view

Bottom view

Stretched grid
• Procedure
1st: stretch grid focused on the north pole. Fig. 1: transgrid
2nd: rotation to the target point: Fig. 2: rotategrid

Fig.１

Fig.２: Center (30N, 120E)