# 4. MITgcm Tutorial Example Experiments¶

The full MITgcm distribution comes with a set of pre-configured numerical experiments. Some of these example experiments are tests of individual parts of the model code, but many are fully fledged numerical simulations. Full tutorials exist for a few of the examples, and are documented in sections Section 4.1 - Section 4.2. The other examples follow the same general structure as the tutorial examples. However, they only include brief instructions in text file README. The examples are located in subdirectories under the directory verification. Each example is briefly described below.

## 4.1. Barotropic Gyre MITgcm Example¶

(in directory verification/tutorial_barotropic_gyre/)

This example experiment demonstrates using the MITgcm to simulate a barotropic, wind-forced, ocean gyre circulation. The experiment is a numerical rendition of the gyre circulation problem described analytically by Stommel in 1948 [Sto48] and Munk in 1950 [Mun50], and numerically in Bryan (1963) [Bry63]. Note this tutorial assumes a basic familiarity with ocean dynamics and geophysical fluid dynamics; readers new to the field may which to consult one of the standard texts on these subjects, such as Vallis (2017) [Val17] or Cushman-Roisin and Beckers (2011) [CRB11].

In this experiment the model is configured to represent a rectangular enclosed box of fluid, \(1200 \times 1200\) km in lateral extent. The fluid depth \(D =\) 5 km. The fluid is forced by a zonal wind stress, \(\tau_x\), that varies sinusoidally in the north-south direction and is constant in time. Topologically the grid is Cartesian and the Coriolis parameter \(f\) is defined according to a mid-latitude beta-plane equation

where \(y\) is the distance along the ‘north-south’ axis of the simulated domain. For this experiment \(f_{0}\) is set to \(10^{-4}\text{s}^{-1}\) and \(\beta = 10^{-11}\text{s}^{-1}\text{m}^{-1}\).

The sinusoidal wind-stress variations are defined according to

where \(L_{y}\) is the lateral domain extent and \(\tau_0\) is set to \(0.1 \text{N m}^{-2}\).

Figure 4.1 summarizes the configuration simulated.

### 4.1.1. Equations Solved¶

The model is configured in hydrostatic form (the MITgcm default). The implicit free surface form of the pressure equation described in Marshall et al. (1997) [MHPA97] is employed. A horizontal Laplacian operator \(\nabla_{h}^2\) provides viscous dissipation. The wind-stress momentum input is added to the momentum equation for the ‘zonal flow’, \(u\). This effectively yields an active set of equations for this configuration as follows:

where \(u\) and \(v\) are the \(x\) and \(y\) components of the flow vector \(\vec{u}\), \(\eta\) is the free surface height, \(A_{h}\) the horizontal Laplacian viscosity, \(\rho_{c}\) is the fluid density, and \(g\) the acceleration due to gravity.

### 4.1.2. Discrete Numerical Configuration¶

The domain is discretized with a uniform grid spacing in the horizontal set to \(\Delta x=\Delta y=20\) km, so that there are sixty grid cells in the \(x\) and \(y\) directions. Vertically the model is configured using a single layer in depth, \(\Delta z\), of 5000 m.

#### 4.1.2.1. Numerical Stability Criteria¶

Let’s start with our choice for the model’s time step. To minimize the amount of required computational resources, typically one opts for as large a time step as possible while keeping the model solution stable. The advective Courant–Friedrichs–Lewy (CFL) condition (see Adcroft 1995 [Adc95]) for an extreme maximum horizontal flow speed is:

The 2 factor on the left is because we have a 2D problem
(in contrast with the more familiar 1D canonical stability analysis); the right hand side is 0.5
due to our default use of Adams-Bashforth2 (see Section 2.5) rather than the more familiar
value of 1 that one would obtain using a forward Euler scheme.
In our configuration, let’s assume our solution will achieve a maximum \(| u | = 1\) ms^{–1}
(in reality, current speeds in our solution will be much smaller). To keep \(\Delta t\) safely
below the stability threshold, let’s choose \(\Delta t\) = 1200 s (= 20 minutes), which
results in \(S_{a}\) = 0.12.

The numerical stability for inertial oscillations using Adams-Bashforth2 (Adcroft 1995 [Adc95])

evaluates to 0.12 for our choice of \(\delta t\), which is below the stability threshold.

There are two general rules in choosing a horizontal Laplacian eddy viscosity \(A_{h}\):

- the resulting Munk layer width should be at least as large (preferably, larger) than the lateral grid spacing;
- the viscosity should be sufficiently small that the model is stable for horizontal friction, given the time step.

Let’s use this first rule to make our choice for \(A_{h}\), and check this value using the second rule. The theoretical Munk boundary layer width (as defined by the solution zero-crossing, see Pedlosky 1987 [Ped87]) is given by:

For our configuration we will choose to resolve a boundary layer of \(\approx\) 100 km,
or roughly across five grid cells, so we set \(A_{h} = 400\) m^{2} s^{–1}
(more precisely, this sets the full width at \(M_{w}\) = 124 km). This choice ensures
that the frictional boundary layer is well resolved.

Given our choice of \(\Delta t\), the stability parameter for the horizontal Laplacian friction (Adcroft 1995 [Adc95])

evaluates to 0.0096, which is well below the stability threshold. As in (4.4) the above criteria is for a 2D problem using Adams-Bashforth2 time stepping, with the 0.6 value on the right replacing the more familiar 1 that is obtained using a forward Euler scheme.

See Section 2.5 for additional details on Adams-Bashforth time-stepping and numerical stability criteria.

### 4.1.3. Code Configuration¶

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

The experiment files

- verification/tutorial_barotropic_gyre/code/SIZE.h
- verification/tutorial_barotropic_gyre/input/data
- verification/tutorial_barotropic_gyre/input/data.pkg
- verification/tutorial_barotropic_gyre/input/eedata
- verification/tutorial_barotropic_gyre/input/bathy.bin
- verification/tutorial_barotropic_gyre/input/windx_cosy.bin

contain the code customizations and parameter settings for this experiment. Below we describe these customizations in detail.

Note: MITgcm’s defaults are configured to simulate an ocean rather than an atmosphere, with vertical \(z\)-coordinates. To model the ocean using pressure coordinates using MITgcm, additional parameter changes are required; see tutorial ocean_in_p. To switch parameters to model an atmosphere, see tutorial Held_Suarez.

#### 4.1.3.1. 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 = 62,
& sNy = 62,
& OLx = 2,
& OLy = 2,
& nSx = 1,
& nSy = 1,
& nPx = 1,
& nPy = 1,
& Nx = sNx*nSx*nPx,
& Ny = sNy*nSy*nPy,
& Nr = 1)
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 )
``` |

Here we show a modified model/inc source code file, customizing MITgcm’s array sizes to our model domain. This file must be uniquely configured for any model setup; using the MITgcm default model/inc/SIZE.h will in fact cause a compilation error. Note that MITgcm’s storage arrays are allocated as static variables (hence their size must be declared in the source code), in contrast to some model codes which declare array sizes dynamically, i.e., through runtime (namelist) parameter settings.

For this first tutorial, our setup and run environment is the most simple possible: we run on a single process (i.e., NOT MPI and NOT multi-threaded) using a single model “tile”. For a more complete explanation of the parameter choices to use multiple tiles, see the tutorial Baroclinic Gyre.

These lines set parameters sNx and sNy, the number of grid points in the \(x\) and \(y\) directions, respectively.

45 46

& sNx = 62, & sNy = 62,

These lines set parameters OLx and OLy in the \(x\) and \(y\) directions, respectively. These values are the overlap extent of a model tile, the purpose of which will be explained in later tutorials. Here, we simply specify the required minimum value (2) in both \(x\) and \(y\).

47 48

& OLx = 2, & OLy = 2,

These lines set parameters nSx, nSy, nPx, and nPy, the number of model tiles and the number of processes in the \(x\) and \(y\) directions, respectively. As discussed above, in this tutorial we configure a single model tile on a single process, so these parameters are all set to the value one.

49 50 51 52

& nSx = 1, & nSy = 1, & nPx = 1, & nPy = 1,

This line sets parameter Nr, the number of points in the vertical dimension. Here we use just a single vertical level.

55

& Nr = 1)

Note these lines summarize the horizontal size of the model domain (

**NOT**to be edited).53 54

& Nx = sNx*nSx*nPx, & Ny = sNy*nSy*nPy,

Further information and examples about how to configure model/inc/SIZE.h are given in Section 6.3.1.

#### 4.1.3.2. 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 | ```
# Model parameters
# Continuous equation parameters
&PARM01
viscAh=4.E2,
f0=1.E-4,
beta=1.E-11,
rhoNil=1000.,
gBaro=9.81,
rigidLid=.FALSE.,
implicitFreeSurface=.TRUE.,
# momAdvection=.FALSE.,
tempStepping=.FALSE.,
saltStepping=.FALSE.,
&
# Elliptic solver parameters
&PARM02
cg2dTargetResidual=1.E-7,
cg2dMaxIters=1000,
&
# Time stepping parameters
&PARM03
nIter0=0,
nTimeSteps=10,
deltaT=1200.0,
pChkptFreq=31104000.0,
chkptFreq=15552000.0,
dumpFreq=15552000.0,
monitorFreq=1200.,
monitorSelect=2,
#-for longer run (3.0 yr):
# nTimeSteps=77760,
# monitorFreq=864000.,
&
# Gridding parameters
&PARM04
usingCartesianGrid=.TRUE.,
delX=62*20.E3,
delY=62*20.E3,
xgOrigin=-20.E3,
ygOrigin=-20.E3,
delR=5000.,
&
# Input datasets
&PARM05
bathyFile='bathy.bin'
zonalWindFile='windx_cosy.bin',
#zonalWindFile='windx_siny.bin',
meridWindFile=,
&
``` |

This file, reproduced completely above, specifies the main parameters for the experiment. The parameters that are significant for this configuration (shown with line numbers to left) are as follows.

##### PARM01 - Continuous equation parameters¶

This line sets parameter viscAh, the horizontal Laplacian viscosity, to \(400\) m

^{2}s^{–1}.4

viscAh=4.E2,

These lines set \(f_0\) and \(\beta\) (the Coriolis parameter f0 and the gradient of the Coriolis parameter beta) for our beta-plane to \(1 \times 10^{-4}\) s

^{–1}and \(1 \times 10^{-11}\) m^{–1}s^{–1}, respectively.5 6

f0=1.E-4, beta=1.E-11,

This line sets parameter rhoNil, a reference density which will also be used as \(\rho_c\) (parameter rhoConst) in (4.1), to 1000 kg/m

^{3}.7

rhoNil=1000.,

This line sets parameter gBaro, the acceleration due to gravity \(g\) (in the free surface terms in (4.1) and (4.2)), to 9.81 m/s

^{2}. This is the MITgcm default value, i.e., the value used if this line were not included in`data`

. One might alter this parameter for a reduced gravity model, or to simulate a different planet, for example.8

gBaro=9.81,

These lines set parameters rigidLid and implicitFreeSurface in order to suppress the rigid lid formulation of the surface pressure inverter and activate the implicit free surface formulation.

9 10

rigidLid=.FALSE., implicitFreeSurface=.TRUE.,

This line sets parameter momAdvection to suppress the (non-linear) momentum of advection terms in the momentum equations. However, note the

`#`

in column 1: this “comments out” the line, so using the above data file verbatim will in fact include the momentum advection terms (i.e., MITgcm default for this parameter is TRUE). We’ll explore the linearized solution (i.e., by removing the leading`#`

) in Section 4.1.5. Note the ability to comment out a line in a namelist file is not part of standard Fortran, but this feature is implemented for all MITgcm namelist files.11

# momAdvection=.FALSE.,

These lines set parameters tempStepping and saltStepping to suppress MITgcm’s forward time integration of temperature and salt in the tracer equations, as these prognostic variables are not relevant for the model solution in this configuration. By default, MITgcm solves equations governing these two (active) tracers; later tutorials will demonstrate how additional passive tracers can be included in the solution. The advantage of NOT solving the temperature and salinity equations is to eliminate many unnecessary computations. In most typical configurations however, one will want the model to compute a solution for \(T\) and \(S\), which typically comprises the majority of MITgcm’s processing time.

12 13

tempStepping=.FALSE., saltStepping=.FALSE.,

##### PARM02 - Elliptic solver parameters¶

The first line sets the tolerance (parameter cg2dTargetResidual) that the 2D conjugate gradient solver, the iterative method used in the pressure method algorithm, will use to test for convergence. The second line sets parameter cg2dMaxIters, the maximum number of iterations. The solver will iterate until the residual falls below this target value (here, set to \(1 \times 10^{-7}\)) or until this maximum number of solver iterations is reached (here, set to a maximum 1000 iterations). Typically, the solver will converge in far fewer than 1000 iterations, but it does not hurt to allow for a large number. The chosen value for the target residual happens to be the MITgcm default, and will serve well in most model configurations.

18 19

cg2dTargetResidual=1.E-7, cg2dMaxIters=1000,

##### PARM03 - Time stepping parameters¶

This line sets the starting (integer) iteration number for the run. Here we set the value to zero, which starts the model from a new, initialized state. If nIter0 is non-zero, the model would require appropriate pickup files (i.e., restart files) in order to continue integration of an existing run.

24

nIter0=0,

This line sets parameter nTimeSteps, the (integer) number of timesteps the model will integrate forward. Below, we have set this to integrate for just 10 time steps, for MITgcm automated testing purposes (Section 5.5). To integrate the solution to near steady state, uncomment the line a few lines further down where we set the value to 77760 time steps. When you make this change, be sure to also comment out the line that sets monitorFreq (see below).

25

nTimeSteps=10,

This line sets parameter deltaT, the timestep used in stepping forward the model, to 1200 seconds. In combination with the larger value of nTimeSteps mentioned above, we have effectively set the model to integrate forward for \(77760 \times 1200 \text{ s} = 3.0\) years (based on 360-day years), long enough for the solution to approach equilibrium.

26

deltaT=1200.0,

These lines control the frequency at which restart (a.k.a. pickup) files are dumped by MITgcm. Here the value of pChkptFreq is set to 31,104,000 seconds (=1.0 years) of model time; this controls the frequency of “permanent” checkpoint pickup files. With permanent files, the model’s iteration number is part of the file name (as the filename “suffix”; see Section 4.1.4.2) in order to save it as a labelled, permanent, pickup state. The value of ChkptFreq is set to 15,552,000 seconds (=0.5 years); the pickup files written at this frequency but will NOT include the iteration number in the filename, instead toggling between

`ckptA`

and`ckptB`

in the filename, and thus these files will be overwritten with new data every 2 \(\times\) 15,552,000 seconds. Temporary checkpoint files can be written more frequently without requiring additional disk space, for example to peruse (or re-run) the model state prior to an instability, or restart following a computer crash, etc. Either type of checkpoint file can be used to restart the model.27 28

pChkptFreq=31104000.0, chkptFreq=15552000.0,

This line sets parameter dumpFreq, frequency of writing model state snapshot diagnostics (of relevance in this setup: variables \(u\), \(v\), and \(\eta\)). Here, we opt for a snapshot of model state every 15,552,000 seconds (=0.5 years), or after every 12960 time steps of integration.

29

dumpFreq=15552000.0,

These lines are set to dump monitor output (see Section 9.4) every 1200 seconds (i.e., every time step) to standard output. While this monitor frequency is needed for MITgcm automated testing, this is too much output for our tutorial run. Comment out this line and uncomment the line where monitorFreq is set to 864,000 seconds, i.e., output every 10 days. Parameter monitorSelect is set to 2 here to reduce output of non-applicable quantities for this simple example.

30 31

monitorFreq=1200., monitorSelect=2,

##### PARM04 - Gridding parameters¶

This line sets parameter usingCartesianGrid, which specifies that the simulation will use a Cartesian coordinate system.

39

usingCartesianGrid=.TRUE.,

These lines set the horizontal grid spacing of the model grid, as vectors delX and delY (i.e., \(\Delta x\) and \(\Delta y\) respectively). This syntax indicates that we specify 62 values in both the \(x\) and \(y\) directions, which matches the domain size as specified in SIZE.h. Grid spacing is set to \(20 \times 10^{3}\) m (=20 km).

40 41

delX=62*20.E3, delY=62*20.E3,

The cartesian grid default origin is (0,0) so here we set the origin with parameters xgOrigin and ygOrigin to (-20000,-20000), accounting for the bordering solid wall. The centers of the grid boxes will thus be at -10 km, 10 km, 30 km, 50 km, …, in both \(x\) and \(y\) directions.

42 43

xgOrigin=-20.E3, ygOrigin=-20.E3,

This line sets parameter delR, the vertical grid spacing in the \(z\)-coordinate (i.e., \(\Delta z\)), to 5000 m.

44

delR=5000.,

##### PARM05 - Input datasets¶

This line sets parameter bathyFile, the name of the bathymetry file. See Section 4.1.3.5 for information about the file format.

49

bathyFile='bathy.bin'

These lines specify the names of the files from which the surface wind stress is read. There is a separate file for the \(x\)-direction (zonalWindFile) and the \(y\)-direction (meridWindFile). Note, here we have left the latter parameter blank, as there is no meridional wind stress forcing in our example.

50 51 52

zonalWindFile='windx_cosy.bin', #zonalWindFile='windx_siny.bin', meridWindFile=,

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

1 2 3 | ```
# Packages
&PACKAGES
&
``` |

This file does not set any namelist parameters, yet is necessary to run – only standard packages (i.e., those compiled in MITgcm by default) are required for this setup, so no other customization is necessary. We will demonstrate how to include additional packages in other tutorial experiments.

#### 4.1.3.4. File input/eedata¶

1 2 3 4 5 6 7 8 9 10 11 | ```
# Example "eedata" file
# Lines beginning "#" are comments
# nTx :: No. threads per process in X
# nTy :: No. threads per process in Y
# debugMode :: print debug msg (sequence of S/R calls)
&EEPARMS
nTx=1,
nTy=1,
&
# Note: Some systems use & as the namelist terminator (as shown here).
# Other systems use a / character.
``` |

This file uses standard default values (i.e., MITgcm default is single-threaded) and does not contain customizations for this experiment.

#### 4.1.3.5. File `input/bathy.bin`

¶

This file is a 2D(\(x,y\)) map of bottom bathymetry, specified as the \(z\)-coordinate of the solid bottom boundary. Here, the value is set to -5000 m everywhere except along the N, S, E, and W edges of the array, where the value is set to 0 (i.e., “land”). The domain in MITgcm is assumed doubly periodic (i.e., periodic in both \(x\)- and \(y\)-directions), so boundary walls are necessary to set up our enclosed box domain. The points are ordered from low to high coordinates in both axes (varying fastest in \(x\)), as a raw binary stream of data that is enumerated in the same way as standard MITgcm 2D horizontal arrays. By default, this file is assumed to contain 32-bit (single precision) binary numbers. The matlab program verification/tutorial_barotropic_gyre/input/gendata.m was used to generate this bathymetry file.

#### 4.1.3.6. File `input/windx_cosy.bin`

¶

Similar to file `input/bathy.bin`

, this file is a 2D(\(x,y\))
map of \(\tau_{x}\) wind stress values, formatted in the same manner.
The units are Nm^{–2}. Although \(\tau_{x}\) is only a function of \(y\) in this experiment,
this file must still define a complete 2D map in order
to be compatible with the standard code for loading forcing fields
in MITgcm. The matlab program verification/tutorial_barotropic_gyre/input/gendata.m
was used to generate this wind stress file. To run the barotropic jet variation of this tutorial example (see Figure 4.4),
you will in fact need to run this
matlab program to generate the file `input/windx_siny.bin`

.

### 4.1.4. Building and running the model¶

To configure and compile the code (following the procedure described in Section 3.5.1):

```
cd build
../../../tools/genmake2 -mods ../code ««-of my_platform_optionFile»»
make depend
make
cd ..
```

To run the model (following the procedure in Section 3.6):

```
cd run
ln -s ../input/* .
ln -s ../build/mitgcmuv .
./mitgcmuv > output.txt
```

#### 4.1.4.1. Standard output¶

Your run’s standard output file should be similar to verification/tutorial_barotropic_gyre/results/output.txt. The standard output is essentially a log file of the model run. The following information is included (in rough order):

- startup information including MITgcm checkpoint release number and other execution environment information, and a list of activated packages (including all default packages, as well as optional packages).
- the text from all
`data.*`

and other critical files (in our example here, eedata, SIZE.h, data, and data.pkg). - information about the grid and bathymetry, including dumps of all grid variables (only if Cartesian or spherical polar coordinates used, as is the case here).
- all runtime parameter choices used by the model, including all model defaults as well as user-specified parameters.
- monitor statistics at regular intervals (as specified by parameter monitorFreq in data. See Section 9.4).
- output from the 2D conjugate gradient solver. More specifically, statistics from the right-hand side of the elliptic equation – for our linear free-surface, see eq. (2.15) – are dumped for every model time step. If the model solution blows up, these statistics will increase to infinity, so one can see exactly when the problem occurred (i.e., to aid in debugging). Additional solver variables, such as number of iterations and residual, are included with the monitor statistics.
- a summary of end-of-run execution information, including user-, wall- and system-time elapsed during execution, and tile communication statistics. These statistics are provided for the overall run, and also broken down by time spent in various subroutines.

Different setups using non-standard packages and/or different parameter choices will include additional or different output as part of the standard output. It is also possible to select more or less output by changing the parameter debugLevel in data; see (missing doc for pkg debug).

`STDERR.0000`

- if errors (or warnings) occurred during the run, helpful warning and/or error message(s) would appear in this file.

#### 4.1.4.2. Other output files¶

In addition to raw binary data files with `.data`

extension, each binary file has a corresponding `.meta`

file. These plain-text files include
information about the array size, precision (i.e., `float32`

or `float64`

), and if relevant, time information and/or
a list of these fields included in the binary file. The `.meta`

files are used by MITgcm utils when binary data are read.

The following output files are generated:

**Grid Data**: see Section 2.11 for definitions and description of the Arakawa C-grid
staggering of model variables.

`XC`

,`YC`

- grid cell center point locations`XG`

,`YG`

- locations of grid cell vertices`RC`

,`RF`

- vertical cell center and cell faces positions`DXC`

,`DYC`

- grid cell center point separations (Figure 2.6 b)`DXG`

,`DYG`

- separation of grid cell vertices (Figure 2.6 a)`DRC`

,`DRF`

- separation of vertical cell centers and faces, respectively`RAC`

,`RAS`

,`RAW`

,`RAZ`

- areas of the grid “tracer cells”, “southern cells”, “western cells” and “vorticity cells”, respectively (Figure 2.6)`hFacC`

,`hFacS`

,`hFacW`

- fractions of the grid cell in the vertical which are “open” as defined in the center and on the southern and western boundaries, respectively. These variables effectively contain the configuration bathymetric (or topographic) information.`Depth`

- bathymetry depths

All these files contain 2D(\(x,y\)) data except `RC`

, `RF`

, `DRC`

, `DRF`

, which are 1D(\(z\)),
and `hFacC`

, `hFacS`

, `hFacW`

, which contain 3D(\(x,y,z\)) data. Units for the grid files depends on one’s choice of model grid;
here, they are all in given in meters (or \(\text{m}^2\) for areas).

All the 2D grid data files contain `.001.001`

in their filename, e.g., `DXC.001.001.data`

– this is the tile number in `.XXX.YYY`

format.
Here, we have just a single tile in both x and y, so both tile numbers are `001`

. Using multiple tiles, the default is that the local tile grid information
would be output separately for each tile (as an example, see the baroclinic gyre tutorial,
which is set up using multiple tiles), producing multiple files for each 2D grid variable.

**State Variable Snapshot Data**:

`Eta.0000000000.001.001.data, Eta.0000000000.001.001.meta`

- this is
a binary data snapshot of model dynamic variable
etaN (the free-surface height) and its meta file, respectively.
Note the tile number is included in the filename, as is the iteration number `0000000000`

, which is simply the time step
(the iteration number here is referred to as the “suffix” in
MITgcm parlance; there are options to change this suffix to something other than iteration number).
In other words, this is a dump of the free-surface height from the initialized state,
iteration 0; if you load up this data file,
you will see it is all zeroes. More interesting is the free-surface
height after some time steps have occurred. Snapshots are written according
to our parameter choice dumpFreq, here set to 15,552,000 seconds, which is every 12960 time steps.
We will examine the model solutions in Section 4.1.5.
The free-surface height is a 2D(\(x,y\)) field.

Snapshot files exist for other prognostic model variables, in particular
filenames starting with `U`

(uVel),
`V`

(uVel), `T`

(theta), and `S`

(salt);
given our setup, these latter two fields
remain uniform in space and time, thus not very interesting until we
explore a baroclinic gyre setup in tutorial_baroclinic_gyre.
These are all 3D(\(x,y,z\)) fields. The format for the file names is similar
to the free-surface height files. Also dumped are snapshots
of diagnosed vertical velocity `W`

(wVel) (note that in non-hydrostatic
simulations, `W`

is a fully prognostic model variable).

**Checkpoint Files**:

The following pickup files are generated:

`pickup.0000025920.001.001.data`

,`pickup.0000025920.001.001.meta`

, etc. - written at frequency set by pChkptFreq`pickup.ckptA.001.001.data`

,`pickup.ckptA.001.001.meta`

,`pickup.ckptB.001.001.data`

,`pickup.ckptB.001.001.meta`

- written at frequency set by ChkptFreq

**Other Model Output Data**: For completeness, here we list the remaining default output files produced by MITgcm (despite being not particularly informative for this simple setup).

`RhoRef.data, RhoRef.meta`

- this is a 1D(\(z\)) array of reference density. Here we have a single level and have not specified an equation of state relation, thus
the file simply contains our prescribed value rhoNil.

`PHrefC.data, PHrefC.meta, PHrefF.data, PHrefF.meta`

- these are 1D(\(z\)) arrays containing reference
hydrostatic “pressure potential” \(\phi = p/\rho_c\) (see Section 1.3.6),
computed at the (vertical grid) cell centers and cell faces, respectively.
In our setup here, `PHrefC`

is simply \(\frac{\rho_c*g*D/2}{\rho_c}\),
i.e., computed at the midpoint of our single vertical cell.

`PH`

, `PHL`

files - these are a 3D(\(x,y,z\)) field of hydrostatic
\(\phi’\) (including free-surface contribution) at cell centers
and a 2D(\(x,y\)) field of ocean bottom \(\phi’\), respectively, as a function of time.
To obtain full \(\phi(t)\) values, `PHrefC`

should be added to `PH`

,
and `PHrefF`

(\(z\) =bottom) should be added to `PHL`

.

### 4.1.5. Model Solution¶

After running the model for 77,760 time steps (3.0 years), the solution is near equilibrium. Given an approximate timescale of one month for barotropic Rossby waves to cross our model domain, one might expect the solution to require several years to achieve an equilibrium state. The model solution of free-surface height \(\eta\) (proportional to streamfunction) at \(t=\) 3.0 years is shown in Figure 4.2. For further details on this solution, particularly examining the effect of the non-linear terms with increasing Reynolds number, the reader is referred to Pedlosky (1987) [Ped87] section 5.11.

Using matlab for example, visualizing output using the utils/matlab/rdmds.m utility to load the
binary data in `Eta.0000077760.001.001.data`

is as simple as:

```
addpath ../../../utils/matlab/
XC=rdmds('XC'); YC=rdmds('YC');
Eta=rdmds('Eta',77760);
contourf(XC/1000,YC/1000,Eta,[-.04:.01:.04]); colorbar;
colormap((flipud(hot))); set(gca,'XLim',[0 1200]); set(gca,'YLim',[0 1200])
```

or using python (you will need to copy utils/python/MITgcmutils/MITgcmutils/mds.py to your run directory before proceeding):

```
import mds
import matplotlib.pyplot as plt
XC = mds.rdmds('XC'); YC = mds.rdmds('YC')
Eta = mds.rdmds('Eta', 77760)
plt.contourf(XC, YC, Eta, np.linspace(-0.02, 0.05,8), cmap='hot_r')
plt.colorbar(); plt.show()
```

Let’s simplify the example by considering the linear problem where we neglect the advection of momentum terms.
In other words, replace \(\frac{Du}{Dt}\) and \(\frac{Dv}{Dt}\) with
\(\frac{\partial u}{\partial t}\) and \(\frac{\partial v}{\partial t}\), respectively, in in (4.1) and (4.2).
To do so, we uncomment (i.e., remove the leading `#`

) in the
line `# momAdvection=.FALSE.,`

in file `data`

and re-run the model. Any existing output files will be overwritten.

For the linearized equations, the Munk layer (equilibrium) analytical solution is given by:

where \(\delta_m= ( \frac { A_{h} }{ \beta } )^{\frac{1}{3}}\). Figure 4.3 displays the MITgcm output after switching off momentum advection vs. the analytical solution to the linearized equations. Success!

Finally, let’s examine one additional simulation where we change the cosine profile of wind stress forcing to a sine profile.
First, run the matlab script verification/tutorial_barotropic_gyre/input/gendata.m to generate the alternate sine
profile wind stress, and place a copy in your run directory. Then,
in file data,
replace the line `zonalWindFile='windx_cosy.bin’,`

with `zonalWindFile='windx_siny.bin’,`

.

The free surface solution given this forcing is shown in Figure 4.4. Two “half gyres” are separated by a strong jet. We’ll look more at the solution to this “barotropic jet” setup in later tutorial examples.

## 4.2. A Rotating Tank in Cylindrical Coordinates¶

(in directory: verification/rotating_tank/)

This example configuration demonstrates using the MITgcm to simulate a laboratory demonstration using a differentially heated rotating annulus of water. The simulation is configured for a laboratory scale on a \(3^{\circ}\times1\mathrm{cm}\) cyclindrical grid with twenty-nine vertical levels of 0.5cm each. This is a typical laboratory setup for illustration principles of GFD, as well as for a laboratory data assimilation project.

example illustration from GFD lab here