# 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 2^{nd} 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 8^{th} 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:

where \(\vec{\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

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:

with:

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:

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{\mathbf{k}}\) is the vertical unity vector, \(\rm{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\):

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).

### 4.7.3.1. Numerical Stability Criteria¶

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

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_{i} < 1\) upper limit for stability. The advective CFL (Adcroft 1995 [Adc95]) for a extreme maximum horizontal flow speed of \(| \vec{u} | = 90 {\rm m/s}\) and the smallest horizontal grid spacing \(\Delta x = 1.1\times10^5 {\rm m}\):

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])

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.

### 4.7.4.1. File input/data¶

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | ```
# ====================
# | Model parameters |
# ====================
#
# Continuous equation parameters
&PARM01
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,
sRef=20*0.0,
no_slip_sides=.FALSE.,
no_slip_bottom=.FALSE.,
buoyancyRelation='ATMOSPHERIC',
eosType='IDEALG',
rotationPeriod=86400.,
implicitFreeSurface=.TRUE.,
exactConserv=.TRUE.,
nonlinFreeSurf=4,
select_rStar=2,
hFacInf=0.2,
hFacSup=2.0,
uniformLin_PhiSurf=.FALSE.,
#hFacMin=0.2,
saltStepping=.FALSE.,
momViscosity=.FALSE.,
vectorInvariantMomentum=.TRUE.,
staggerTimeStep=.TRUE.,
readBinaryPrec=64,
writeBinaryPrec=64,
&
# Elliptic solver parameters
&PARM02
cg2dMaxIters=200,
#cg2dTargetResidual=1.E-12,
cg2dTargetResWunit=1.E-17,
&
# Time stepping parameters
&PARM03
deltaT=450.,
#nIter0=276480,
startTime=124416000.,
#- run for 1 year (192.iterations x 450.s = 1.day, 360*192=69120):
#nTimeSteps=69120,
#forcing_In_AB=.FALSE.,
tracForcingOutAB=1,
abEps=0.1,
pChkptFreq=31104000.,
chkptFreq=2592000.,
dumpFreq=2592000.,
#monitorFreq=43200.,
taveFreq=0.,
#- to run a short test (2.h):
nTimeSteps=16,
monitorFreq=1.,
&
# Gridding parameters
&PARM04
usingCurvilinearGrid=.TRUE.,
horizGridFile='grid_cs32',
radius_fromHorizGrid=6370.E3,
delR=20*50.E2,
&
# Input datasets
&PARM05
#topoFile='topo.cs.bin',
&
``` |

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,

no_slip_sides=.FALSE.,

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,

no_slip_bottom=.FALSE.,

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,

buoyancyRelation='ATMOSPHERIC',

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,

eosType='IDEALG',

Selects the ideal gas equation of state.

Line 15,

implicitFreeSurface=.TRUE.,

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

Line 16,

exactConserv=.TRUE.,

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,

nonlinFreeSurf=4, select_rStar=2,

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,

uniformLin_PhiSurf=.FALSE.,

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_{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_{Ref}\).Line 23-24,

saltStepping=.FALSE., momViscosity=.FALSE.,

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

Line 25,

vectorInvariantMomentum=.TRUE.,

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

Line 26,

staggerTimeStep=.TRUE.,

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

Lines 27-28,

readBinaryPrec=64, writeBinaryPrec=64,

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

Line 33,

cg2dMaxIters=200,

Sets maximum number of iterations the 2-D conjugate gradient solver will use,

**irrespective of convergence criteria being met**.Line 35,

cg2dTargetResWunit=1.E-17,

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,

deltaT=450.,

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

Line 42,

startTime=124416000.,

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,

#nTimeSteps=69120,

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

Lines 54-55,

nTimeSteps=16, monitorFreq=1.,

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,

pChkptFreq=31104000.,

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,

chkptFreq=2592000.,

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,

dumpFreq=2592000.,

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

Line 51,

#monitorFreq=43200.,

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,

usingCurvilinearGrid=.TRUE.,

Set the horizontal type of grid to curvilinear grid.

Line 61,

horizGridFile='grid_cs32',

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,

delR=20*50.E2,

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,

#topoFile='topo.cs.bin'

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..

### 4.7.4.2. File input/data.pkg¶

1 2 3 4 5 6 | ```
# Packages
&PACKAGES
useSHAP_FILT=.TRUE.,
useDiagnostics=.TRUE.,
#useMNC=.TRUE.,
&
``` |

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,

useSHAP_FILT=.TRUE.,

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

Line 4,

useDiagnostics=.TRUE.,

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

Line 5,

#useMNC=.TRUE.,

This line would select pkg/mnc for I/O but is commented out.

### 4.7.4.3. File input/data.shap¶

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | ```
# Shapiro Filter parameters
&SHAP_PARM01
shap_filt_uvStar=.FALSE.,
shap_filt_TrStagg=.TRUE.,
Shap_funct=2,
nShapT=0,
nShapUV=4,
#nShapTrPhys=0,
nShapUVPhys=4,
#Shap_TrLength=140000.,
#Shap_uvLength=110000.,
#Shap_Trtau=5400.,
#Shap_uvtau=1800.,
#Shap_diagFreq=2592000.,
&
``` |

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,

Shap_funct=2,

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

Lines 6-7,

nShapT=0, nShapUV=4,

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,

nShapUVPhys=4,

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,

#Shap_Trtau=5400., #Shap_uvtau=1800.,

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.

### 4.7.4.4. File input/eedata¶

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

This file uses standard default values except line 6:

```
useCubedSphereExchange=.TRUE.,
```

This line selects the cubed-sphere specific exchanges to to connect tiles and faces edges.

### 4.7.4.5. File code/SIZE.h¶

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

Four lines are customized in this file for the current experiment

Line 45,

sNx=32,

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

Line 46,

sNy=32,

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

Line 49,

nSx=6,

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,

Nr=20,

sets the vertical domain extent in grid points.

### 4.7.4.6. File code/packages.conf¶

1 2 3 4 5 6 | ```
#-- list of packages (or group of packages) to compile for this experiment:
exch2
gfd
shap_filt
diagnostics
mnc
``` |

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,

gfd

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

Line 2,

shap_filt

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

Line 3,

exch2

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,

diagnostics

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

Line 5,

mnc

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

### 4.7.4.7. File code/CPP_OPTIONS.h¶

This file uses the standard default except for:

```
#define NONLIN_FRSURF
```

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

### 4.7.4.8. Other Files¶

Other files relevant to this experiment are

`input/grid_cs32.face00[n].bin`

, with \(n=1,2,3,4,5,6\)

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)