4.7. Held-Suarez Atmosphere

(in directory: verification/tutorial_held_suarez_cs/)

This example illustrates the use of the MITgcm as an atmospheric GCM, using simple Held and Suarez (1994) [HS94] forcing to simulate atmospheric dynamics on global scale. The set-up uses the rescaled pressure coordinate (\(p^*\)) of Adcroft and Campin (2004) [AC04] in the vertical direction, with 20 equally-spaced levels, and the conformal cube-sphere grid (C32) described in Adcroft et al. (2004) [ACHM04].

4.7.1. Overview

This example demonstrates using the MITgcm to simulate the planetary atmospheric circulation, with flat orography and simplified forcing. In particular, only dry air processes are considered and radiation effects are represented by a simple Newtonian cooling, Thus, this example does not rely on any particular atmospheric physics package. This kind of simplified atmospheric simulation has been widely used in GFD-type experiments and in intercomparison projects of AGCM dynamical cores (Held and Suarez 1994 [HS94]).

The horizontal grid is obtain from the projection of a uniform gridded cube to the sphere. Each of the 6 faces has the same resolution, with 32 \(\times\) 32 grid points. The equator coincides with a grid line and crosses through the middle in 4 of the 6 faces, leaving 2 faces for the northern and southern polar regions. This curvilinear grid requires the use of the 2nd generation exchange topology (pkg/exch2) to connect tile and face edges, but without any limitation on the number of processors.

The use of the \(p^*\) coordinate with 20 equally spaced levels (20 \(\times\) 50 mb, from \(p^*=1000\) mb to \(0\) at the top of the atmosphere) follows the choice of Held and Suarez (1994) [HS94]. Note that without topography, the \(p^*\) coordinate and the normalized pressure coordinate (\(\sigma_p\)) coincide exactly. Both viscosity and diffusivity are set to zero here, but an 8th order Shapiro (1970) [Sha70] filter is applied to both momentum and potential temperature, to remove selectively grid scale noise. Apart from the horizontal grid, this experiment is made very similar to the grid-point model case used in the Held and Suarez (1994) [HS94] study.

At this resolution, the configuration can be integrated forward for many years on a single processor desktop computer.

4.7.2. Forcing

The model is forced by relaxation to a radiative equilibrium temperature from Held and Suarez (1994) [HS94]. A linear frictional drag (Rayleigh damping) is applied in the lower part of the atmosphere and accounts for surface friction and momentum dissipation in the boundary layer. Altogether, this yields the following forcing from Held and Suarez (1994) [HS94] that is applied to the fluid:

(4.48)\[\vec{\boldsymbol{\cal F}_\mathbf{v}} = -k_\mathbf{v}(p)\vec{\mathbf{v}}_h\]
(4.49)\[{\cal F}_{\theta} = -k_{\theta}(\varphi,p)[\theta-\theta_{eq}(\varphi,p)]\]

where \(\vec{\boldsymbol{\cal F}_\mathbf{v}}\), \({\cal F}_{\theta}\), are the forcing terms in the zonal and meridional momentum and in the potential temperature equations, respectively. The term \(k_\mathbf{v}\) in (4.48) applies a Rayleigh damping that is active within the planetary boundary layer. It is defined so as to decay as pressure decreases according to

\[\begin{split}\begin{aligned} \label{eq:eg-hs-define_kv} k_\mathbf{v} & = k_{f}~\max[0,~(p^*/P^{0}_{s}-\sigma_{b})/(1-\sigma_{b})] \\ \sigma_{b} & = 0.7 ~~{\rm and}~~ k_{f} = 1/86400 ~{\rm s}^{-1}\end{aligned}\end{split}\]

where \(p^*\) is the pressure level of the cell center and \(P^{0}_{s}\) is the pressure at the base of the atmospheric column, which is constant and uniform here (\(= 10^5 {\rm Pa}\)), in the absence of topography.

The equilibrium temperature \(\theta_{eq}\) and relaxation time scale \(k_{\theta}\) are set to:

(4.50)\[\theta_{eq}(\varphi,p^*) = \max \{ 200 (P^{0}_{s}/p^*)^\kappa, \nonumber \hspace{2mm} 315 - \Delta T_y~\sin^2(\varphi) - \Delta \theta_z \cos^2(\varphi) \log(p^*/P^{0}_s) \}\]
(4.51)\[k_{\theta}(\varphi,p^*) = k_a + (k_s -k_a)~\cos^4(\varphi)~\max \{ 0,\nonumber \hspace{2mm} (p^*/P^{0}_{s}-\sigma_{b})/(1-\sigma_{b}) \}\]


\[\begin{split}\begin{aligned} \Delta T_y = 60 \text{ K, } k_a &= 1/(40 \cdot 86400) ~{\rm s}^{-1}\\ \Delta \theta_z = 10 \text{ K, } k_s &= 1/(4 \cdot 86400) ~{\rm s}^{-1}\end{aligned}\end{split}\]

Initial conditions correspond to a resting state with horizontally uniform stratified fluid. The initial temperature profile is simply the horizontal average of the radiative equilibrium temperature.

4.7.3. Set-up description

The model is configured in hydrostatic form, using non-Boussinesq \(p^*\) coordinate. The vertical resolution is uniform, \(\Delta p^* = 50 \times 10^2\) Pa, with 20 levels, from \(p^*=10^5\) Pa to \(0\) at the top. The domain is discretized using the C32 cube-sphere grid (see Adcroft et al. 2004 [ACHM04]) that covers the whole sphere with a relatively uniform grid spacing. The resolution at the equator or along the Greenwich meridian is similar to a 128 \(\times\) 64 equally spaced longitude-latitude grid, but requires 25% less grid points. Grid spacing and grid-point location are not computed by the model, but instead read from files.

The vector-invariant form of the momentum equation (see Section 2.15) is used so that no explicit metrics are necessary.

Applying the vector-invariant discretization to the atmospheric equations (1.59), and adding the forcing terms (4.48), (4.49) on the right-hand-side, leads to the set of equations that are solved in this configuration:

(4.52)\[\frac{\partial \vec{\mathbf{v}}_h}{\partial t} +(f + \zeta)\hat{\boldsymbol{k}} \times \vec{\mathbf{v}}_h + \nabla_{p} (\mathrm{KE}) + \omega \frac{\partial \vec{\mathbf{v}}_h }{\partial p} + \nabla_p \Phi ^{\prime } = -k_\mathbf{v} \vec{\mathbf{v}}_h\]
\[\frac{\partial \Phi ^{\prime }}{\partial p} +\frac{\partial \Pi }{\partial p}\theta ^{\prime } =0\]
\[\nabla _{p}\cdot \vec{\mathbf{v}}_h+\frac{\partial \omega }{ \partial p} =0\]
\[\frac{\partial \theta }{\partial t} + \nabla _{p}\cdot (\theta \vec{\mathbf{v}}_h) + \frac{\partial (\theta \omega)}{\partial p} = -k_{\theta}[\theta-\theta_{\rm eq}]\]

where \(\vec{\mathbf{v}}_h\) and \(\omega = \frac{Dp}{Dt}\) are the horizontal velocity vector and the vertical velocity in pressure coordinate, \(\zeta\) is the relative vorticity and \(f\) the Coriolis parameter, \(\hat{\boldsymbol{k}}\) is the vertical unity vector, \(\mathrm{KE}\) is the kinetic energy, \(\Phi\) is the geopotential, and \(\Pi\) the Exner function (\(\Pi = C_p (p/p_c)^\kappa\) with \(p_c = 10^5\) Pa). Primed variables correspond to anomaly from the resting, uniformly stratified state.

As described in Section 2, the continuity equation is integrated vertically, to give a prognostic equation for the surface pressure \(p_s\):

\[\frac{\partial p_s}{\partial t} + \nabla_{h}\cdot \int_{0}^{p_s} \vec{\mathbf{v}}_h dp = 0\]

The implicit free surface form of the pressure equation described in Marshall et al. (1997) [MHPA97] is employed to solve for \(p_s\); Vertically integrating the hydrostatic balance gives the geopotential \(\Phi'\) and allows one to step forward the momentum equation (4.52). The potential temperature, \(\theta\), is stepped forward using the new velocity field (see Section 2.8). Numerical Stability Criteria

The numerical stability for inertial oscillations (see Adcroft 1995 [Adc95]):

(4.53)\[S_{\rm inert} = f^{2} {\Delta t}^2\]

evaluates to \(4 \times10^{-3}\) at the poles, for \(f=2\Omega\sin(\pi / 2) = 1.45\times10^{-4}~{\rm s}^{-1}\), which is well below the \(S_{\rm inert} < 1\) upper limit for stability. The advective CFL (Adcroft 1995 [Adc95]) for a extreme maximum horizontal flow speed of \(| \vec{\mathbf{u}} | = 90 {\rm m/s}\)  and the smallest horizontal grid spacing \(\Delta x = 1.1\times10^5 {\rm m}\):

(4.54)\[S_{\rm adv} = \frac{| \vec{\mathbf{u}} | \Delta t}{ \Delta x}\]

evaluates to 0.37, which is close to the stability limit of 0.5. The stability parameter for internal gravity waves propagating with a maximum speed of \(c_{g}=100~{\rm m/s}\) (Adcroft 1995 [Adc95])

(4.55)\[S_{c} = \frac{c_{g} \Delta t}{ \Delta x}\]

evaluates to \(4 \times 10^{-1}\). This is close to the linear stability limit of 0.5.

4.7.4. Experiment Configuration

The model configuration for this experiment resides under the directory verification/tutorial_held_suarez_cs/. The experiment files

contain the code customizations and parameter settings for these experiments. Below we describe the customizations to these files associated with this experiment. File input/data

Listing 4.26 verification/tutorial_held_suarez_cs/input/data
 1# ====================
 2# | Model parameters |
 3# ====================
 5# Continuous equation parameters
 6 &PARM01
 7 tRef=295.2, 295.5, 295.9, 296.3, 296.7, 297.1, 297.6, 298.1, 298.7, 299.3,
 8      300.0, 300.7, 301.9, 304.1, 308.0, 315.1, 329.5, 362.3, 419.2, 573.8,
 9 sRef=20*0.0,
10 no_slip_sides=.FALSE.,
11 no_slip_bottom=.FALSE.,
12 buoyancyRelation='ATMOSPHERIC',
13 eosType='IDEALG',
14 rotationPeriod=86400.,
15 implicitFreeSurface=.TRUE.,
16 exactConserv=.TRUE.,
17 nonlinFreeSurf=4,
18 select_rStar=2,
19 hFacInf=0.2,
20 hFacSup=2.0,
21 uniformLin_PhiSurf=.FALSE.,
23 saltStepping=.FALSE.,
24 momViscosity=.FALSE.,
25 vectorInvariantMomentum=.TRUE.,
26 staggerTimeStep=.TRUE.,
27 readBinaryPrec=64,
28 writeBinaryPrec=64,
29 &
31# Elliptic solver parameters
32 &PARM02
33 cg2dMaxIters=200,
35 cg2dTargetResWunit=1.E-17,
36 &
38# Time stepping parameters
39 &PARM03
40 deltaT=450.,
42 startTime=124416000.,
43#- run for 1 year (192.iterations x 450.s = 1.day, 360*192=69120):
46 tracForcingOutAB=1,
47 abEps=0.1,
48 pChkptFreq=31104000.,
49 chkptFreq=2592000.,
50 dumpFreq=2592000.,
52 taveFreq=0.,
53#- to run a short test (2.h):
54 nTimeSteps=16,
55 monitorFreq=1.,
56 &
58# Gridding parameters
59 &PARM04
60 usingCurvilinearGrid=.TRUE.,
61 horizGridFile='grid_cs32',
62 radius_fromHorizGrid=6370.E3,
63 delR=20*50.E2,
64 &
66# Input datasets
67 &PARM05
69 &

This file specifies the main parameters for the experiment. The parameters that are significant for this configuration are:

  • Lines 7-8,

    tRef=295.2, 295.5, 295.9, 296.3, 296.7, 297.1, 297.6, 298.1, 298.7, 299.3,
         300.0, 300.7, 301.9, 304.1, 308.0, 315.1, 329.5, 362.3, 419.2, 573.8,

    set reference values for potential temperature (in kelvins) at each model level. The entries are ordered like model level, from surface up to the top. Density is calculated from anomalies at each level evaluated with respect to the reference values set here.

  • Line 10,


    this line selects a free-slip lateral boundary condition for the horizontal Laplacian friction operator, e.g., \(\frac{\partial u}{\partial y}\)=0 along boundaries in \(y\) and \(\frac{\partial v}{\partial x}\)=0 along boundaries in \(x\).

  • Line 11,


    this line selects a free-slip boundary condition at the top, in the vertical Laplacian friction operator, e.g., \(\frac{\partial u}{\partial p} = \frac{\partial v}{\partial p} = 0\).

  • Line 12,


    this line sets the type of fluid and the type of vertical coordinate to use, which, in this case, is air with a pressure-like coordinate (\(p\) or \(p^*\)).

  • Line 13,


    Selects the ideal gas equation of state.

  • Line 15,


    Selects the way the barotropic equation is solved, here using the implicit free-surface formulation.

  • Line 16,


    Explicitly calculate (again) the surface pressure changes from the divergence of the vertically integrated horizontal flow, after the implicit free surface solver and filters are applied.

  • Lines 17-18,


    Select the non-linear free surface formulation, using \(r^*\) vertical coordinate (here \(p^*\)). Note that, except for the default (= 0), other values of those two parameters are only permitted for testing/debugging purpose.

  • Line 21,


    Select the linear relation between surface geopotential anomaly and surface pressure anomaly to be evaluated from \(\frac{\partial \Phi_s}{\partial p_s} = 1/\rho(\theta_{\rm ref})\) (see Section 2.10.2). Note that using the default (=.TRUE.), the constant \(1/\rho_0\) is used instead, and is not necessarily consistent with other parts of the geopotential that rely on \(\theta_{\rm ref}\).

  • Line 23-24,


    Do not step forward water vapor and do not compute viscous terms. This saves computer time.

  • Line 25,


    Select the vector-invariant form to solve the momentum equation.

  • Line 26,


    Select the staggered time-stepping (rather than synchronous time stepping).

  • Lines 27-28,


    Sets format for reading binary input datasets and writing output fields to use 64-bit representation for floating-point numbers.

  • Line 33,


    Sets maximum number of iterations the 2-D conjugate gradient solver will use, irrespective of convergence criteria being met.

  • Line 35,


    Sets the tolerance (in units of \(\omega\)) which the 2-D conjugate gradient solver will use to test for convergence in equation (2.15) to \(1 \times 10^{-17}\) Pa/s. Solver will iterate until tolerance falls below this value or until the maximum number of solver iterations is reached.

  • Line 40,


    Sets the timestep \(\Delta t\) used in the model to 450 seconds (= 1/8 hour).

  • Line 42,


    Sets the starting time, in seconds, for the model time counter. A non-zero starting time requires the initial state read from a pickup file. By default the pickup file is named according to the integer number (nIter0) of time steps in the startTime value (nIter0 = startTime / deltaT).

  • Line 44,


    A commented out setting for the length of the simulation (in number of timesteps) that corresponds to 1-year simulation.

  • Lines 54-55,


    Sets the length of the simulation (in number of timesteps) and the frequency (in seconds) for “monitor” output to 16 iterations and 1 seconds respectively. This choice corresponds to a short simulation test.

  • Line 48,


    Sets the time interval, in seconds, between 2 consecutive “permanent” pickups (“permanent checkpoint frequency”) that are used to restart the simulation, to 1 year.

  • Line 48,


    Sets the time interval, in seconds, between two consecutive “temporary” pickups (“checkpoint frequency”) to one month. The “temporary” pickup file name is alternatively “ckptA” and “ckptB”; these pickups (as opposed to the permanent ones) are designed to be over-written by the model as the simulation progresses.

  • Line 50,


    Set the frequency (in seconds) for the snapshot output to 1 month.

  • Line 51,


    A commented out line setting the frequency (in seconds) for the “monitor” output to 12 h. This frequency fits better with the longer simulation of one year.

  • Line 60,


    Set the horizontal type of grid to curvilinear grid.

  • Line 61,


    Set the root for the grid file name to grid_cs32. The grid-file names are derived from the root, adding a suffix with the face number (e.g., .face001.bin, .face002.bin \(\cdots\) )

  • Lines 63,


    This line sets the increments in pressure units to 20 equally thick levels of \(50 \times 10^2\) Pa each. This defines the origin (interface \(k=1\)) of the vertical pressure axis, with decreasing pressure as the level index \(k\) increases.

  • Line 68,


    This commented out line would set the file name of a 2-D orography file, in units of meters, to topo.cs.bin.

Other lines in the file input/data are standard values that are described in Section 3.. File input/data.pkg

Listing 4.27 verification/tutorial_held_suarez_cs/input/data.pkg
1# Packages
4 useDiagnostics=.TRUE.,
6 &

This file specifies the additional packages that the model uses for the experiment. Note that some packages are used by default (e.g., pkg/generic_advdiff) and some others are selected according to parameter in input/data.pkg (e.g., pkg/mom_vecinv). The additional packages that are used for this configuration are

  • Line 3,


    This line selects the Shapiro filter (Shapiro 1970 [Sha70]) (pkg/shap_filt) to be used in this experiment.

  • Line 4,


    This line selects pkg/diagnostics to be used in this experiment.

  • Line 5,


    This line would select pkg/mnc for I/O but is commented out. File input/data.shap

Listing 4.28 verification/tutorial_held_suarez_cs/input/data.shap
 1# Shapiro Filter parameters
 3 shap_filt_uvStar=.FALSE.,
 4 shap_filt_TrStagg=.TRUE.,
 5 Shap_funct=2,
 6 nShapT=0,
 7 nShapUV=4,
 9 nShapUVPhys=4,
15 &

This file specifies the parameters that the model uses for the Shapiro filter package (Shapiro 1970 [Sha70]), see Section 2.18. The parameters that are significant for this configuration are:

  • Line 5,


    This line selects which Shapiro filter function to use, here S2, for this experiment (see Section 2.18).

  • Lines 6-7,


    These lines select the order of the Shapiro filter for active tracers (\(\theta\) and \(q\)) and momentum (\(u,v\)) respectively. In this case, no filter is applied to active tracers. Regarding the momentum, this sets the integer parameter \(n\) to 4, in the equations of Section 2.18, which corresponds to a 8th order filter.

  • Line 9,


    This line selects the order of the physical space filter (filter function S2g, see Section 2.18) that applies to \(u,v\). The difference nShapUV - nShapUVPhys corresponds to the order of the computational filter (filter function S2c, see Section 2.18).

  • Lines 12-13,


    These commented lines would have set the time scale of the filter (in seconds), for \(\theta\), \(q\) and for \(u\), \(v\) respectively, to 5400 s (90 min) and 1800 s (30 min). Without explicitly setting those timescales, the default is used, which corresponds to the model timestep. File input/eedata

Listing 4.29 verification/tutorial_held_suarez_cs/input/eedata
 1# Example "eedata" file
 2# Lines beginning "#" are comments
 3# nTx - No. threads per process in X
 4# nTy - No. threads per process in Y
 6 useCubedSphereExchange=.TRUE.,
 7# Activate one line below to support 2, 3 or 6 way multi-threading
11 &
12# Note: Some systems use & as the
13# namelist terminator. Other systems
14# use a / character (as shown here).

This file uses standard default values except line 6:


This line selects the cubed-sphere specific exchanges to to connect tiles and faces edges. File code/SIZE.h

Listing 4.30 verification/tutorial_held_suarez_cs/code/SIZE.h
 4C    include SIZE.h
 5C    !DESCRIPTION: \bv
 6C     *==========================================================*
 7C     | SIZE.h Declare size of underlying computational grid.
 8C     *==========================================================*
 9C     | The design here supports a three-dimensional model grid
10C     | with indices I,J and K. The three-dimensional domain
11C     | is comprised of nPx*nSx blocks (or tiles) of size sNx
12C     | along the first (left-most index) axis, nPy*nSy blocks
13C     | of size sNy along the second axis and one block of size
14C     | Nr along the vertical (third) axis.
15C     | Blocks/tiles have overlap regions of size OLx and OLy
16C     | along the dimensions that are subdivided.
17C     *==========================================================*
18C     \ev
20C     Voodoo numbers controlling data layout:
21C     sNx :: Number of X points in tile.
22C     sNy :: Number of Y points in tile.
23C     OLx :: Tile overlap extent in X.
24C     OLy :: Tile overlap extent in Y.
25C     nSx :: Number of tiles per process in X.
26C     nSy :: Number of tiles per process in Y.
27C     nPx :: Number of processes to use in X.
28C     nPy :: Number of processes to use in Y.
29C     Nx  :: Number of points in X for the full domain.
30C     Ny  :: Number of points in Y for the full domain.
31C     Nr  :: Number of points in vertical direction.
33      INTEGER sNx
34      INTEGER sNy
35      INTEGER OLx
36      INTEGER OLy
37      INTEGER nSx
38      INTEGER nSy
39      INTEGER nPx
40      INTEGER nPy
41      INTEGER Nx
42      INTEGER Ny
43      INTEGER Nr
44      PARAMETER (
45     &           sNx =  32,
46     &           sNy =  32,
47     &           OLx =   2,
48     &           OLy =   2,
49     &           nSx =   6,
50     &           nSy =   1,
51     &           nPx =   1,
52     &           nPy =   1,
53     &           Nx  = sNx*nSx*nPx,
54     &           Ny  = sNy*nSy*nPy,
55     &           Nr  =  20)
57C     MAX_OLX :: Set to the maximum overlap region size of any array
58C     MAX_OLY    that will be exchanged. Controls the sizing of exch
59C                routine buffers.
62      PARAMETER ( MAX_OLX = OLx,
63     &            MAX_OLY = OLy )

Four lines are customized in this file for the current experiment

  • Line 45,


    sets the lateral domain extent in grid points along the \(x\)-direction, for one face.

  • Line 46,


    sets the lateral domain extent in grid points along the \(y\)-direction, for one face.

  • Line 49,


    sets the number of tiles in the \(x\)-direction, for the model domain decomposition. In this simple case (single processor, with one tile per face), this number corresponds to the total number of faces.

  • Line 55,


    sets the vertical domain extent in grid points. File code/packages.conf

Listing 4.31 verification/tutorial_held_suarez_cs/input/code/packages.conf
1#-- list of packages (or group of packages) to compile for this experiment:

This file specifies the packages that are compiled and made available for this experiment. The additional packages that are used for this configuration are

  • Line 1,


    This line selects the standard set of packages that are used by default.

  • Line 2,


    This line makes the Shapiro filter package available for this experiment.

  • Line 3,


    This line selects pkg/exch2 to be compiled and used in this experiment. Note that at present, no such parameter useEXCH2 exists and therefore this package is always used when it is compiled.

  • Line 4,


    This line selects pkg/diagnostics to be compiled, and makes it available for this experiment.

  • Line 5,


    This line selects the pkg/mnc to be compiled, and makes it available for this experiment. File code/CPP_OPTIONS.h

This file uses the standard default except for:


This line enables the non-linear free-surface part of the code, which is required for the \(p^*\) coordinate formulation. Other Files

Other files relevant to this experiment are

contain the code customizations and binary input files for this experiment. The file apply_forcing.F contains four subroutines that calculate the forcing terms (i.e., right-hand side terms) in the momentum equation (4.48), EXTERNAL_FORCING_U and EXTERNAL_FORCING_V and in the potential temperature equation (4.49), EXTERNAL_FORCING_T. The water-vapor forcing subroutine (EXTERNAL_FORCING_S) is left empty for this experiment. The grid-files input/grid_cs32.face00[n].bin, with \(n=1,2,3,4,5,6\), are binary files (direct-access, big-endian 64 bit reals) that contains all the cubed-sphere grid lengths, areas and grid-point positions, with one file per face. Each file contains 18 2-D arrays (dimension 33 \(\times\) 33) that correspond to the model variables: XC YC DXF DYF RA XG YG DXV DYU RAZ DXC DYC RAW RAS DXG DYG AngleCS AngleSN (see model/inc/GRID.h)