(in directory :filelink:`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 :cite:`stommel:48` and Munk in 1950 :cite:`munk:50`, and numerically in Bryan (1963) :cite:`bryan:63`. 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) :cite:`vallis:17` or Cushman-Roisin and Beckers (2011) :cite:`cushmanroisin:11`.
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 H = 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
f(y) = f_{0}+\beta y
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
\tau_x(y) = -\tau_{0} \cos\left(\pi \frac{y}{L_y}\right)
where L_{y} is the lateral domain extent and \tau_0 is set to 0.1 \text{N m}^{-2}.
:numref:`baro_simulation_config` summarizes the configuration simulated.
Schematic of simulation domain and wind-stress forcing function for barotropic gyre numerical experiment. The domain is enclosed by solid walls at x= 0, 1200 km and at y= 0, 1200 km.
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) :cite:`marshall:97a` 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:
\frac{Du}{Dt} - fv + g\frac{\partial \eta}{\partial x} - A_{h}\nabla_{h}^2u
= \frac{\tau_{x}}{\rho_{c}H}
\frac{Dv}{Dt} + fu + g\frac{\partial \eta}{\partial y} - A_{h}\nabla_{h}^2v
= 0
\frac{\partial \eta}{\partial t} + \nabla _{h}\cdot (H \vec{\bf u})
= 0
where u and v are the x and y components of the flow vector \vec{\bf 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.
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 ocean grid cells in the x and y directions. The numerical domain includes a border row of "land" cell surrounding the ocean cells, so the numerical grid size is 62\times62 (if these land cells were not included, the domain would be periodic in both the x and y directions).
Vertically the model is configured using a single layer in depth, \Delta z, of 5000 m.
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 :cite:`adcroft:95`) for an extreme maximum horizontal flow speed is:
S_{\rm adv} = 2 \left( \frac{ |u| \Delta t}{ \Delta x} \right) < 0.5 \text{ for stability}
The 2 factor on the left is because we have a 2-D problem (in contrast with the more familiar 1-D canonical stability analysis); the right hand side is 0.5 due to our default use of Adams-Bashforth2 (see :numref:`adams-bashforth`) 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_{\rm adv} = 0.12.
The numerical stability for inertial oscillations using Adams-Bashforth II (Adcroft 1995 :cite:`adcroft:95`)
S_{\rm inert} = f {\Delta t} < 0.5 \text{ for stability}
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 :cite:`pedlosky:87`) is given by:
M = \frac{2\pi}{\sqrt{3}} \left( \frac { A_{h} }{ \beta } \right) ^{\frac{1}{3}}
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 m2 s--1 (more precisely, this sets the full width at M = 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 :cite:`adcroft:95`)
S_{\rm Lh} = 2 \left( 4 \frac{A_{h} \Delta t}{{\Delta x}^2} \right) < 0.6 \text{ for stability}
evaluates to 0.0096, which is well below the stability threshold. As in :eq:`eq_baro_cfl_stability` 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 :numref:`adams-bashforth` for additional details on Adams-Bashforth time-stepping and numerical stability criteria.
The model configuration for this experiment resides under the directory :filelink:`verification/tutorial_barotropic_gyre/`.
The experiment files
- :filelink:`verification/tutorial_barotropic_gyre/code/SIZE.h`
- :filelink:`verification/tutorial_barotropic_gyre/input/data`
- :filelink:`verification/tutorial_barotropic_gyre/input/data.pkg`
- :filelink:`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.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/code/SIZE.h
:linenos:
:caption: verification/tutorial_barotropic_gyre/code/SIZE.h
Here we show a modified :filelink:`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 :filelink:`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 :ref:`"tile" <tile_description>`. For a more complete explanation of the parameter choices to use multiple tiles, see the tutorial Baroclinic Gyre.
These lines set parameters :varlink:`sNx` and :varlink:`sNy`, the number of grid points in the x and y directions, respectively.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/code/SIZE.h :start-at: sNx = :end-at: sNy = :lineno-match:These lines set parameters :varlink:`OLx` and :varlink:`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.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/code/SIZE.h :start-at: OLx = :end-at: OLy = :lineno-match:These lines set parameters :varlink:`nSx`, :varlink:`nSy`, :varlink:`nPx`, and :varlink:`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.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/code/SIZE.h :start-at: nSx = :end-at: nPy = :lineno-match:This line sets parameter :varlink:`Nr`, the number of points in the vertical dimension. Here we use just a single vertical level.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/code/SIZE.h :start-at: Nr = :end-at: Nr = :lineno-match:Note these lines summarize the horizontal size of the model domain (NOT to be edited).
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/code/SIZE.h :start-at: sNx*nSx*nPx :end-at: sNy*nSy*nPy :lineno-match:
Further information and examples about how to configure :filelink:`model/inc/SIZE.h` are given in :numref:`specify_decomp`.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data
:linenos:
:caption: verification/tutorial_barotropic_gyre/input/data
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.
This line sets parameter :varlink:`viscAh`, the horizontal Laplacian viscosity, to 400 m2 s--1.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: viscAh :end-at: viscAh :lineno-match:These lines set f_0 and \beta (the Coriolis parameter :varlink:`f0` and the gradient of the Coriolis parameter :varlink:`beta`) for our beta-plane to 1 \times 10^{-4} s--1 and 1 \times 10^{-11} m--1s--1, respectively.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: f0 :end-at: beta :lineno-match:This line sets parameter :varlink:`rhoConst`, the Boussinesq reference density \rho_c in :eq:`baro_model_eq_u`, to 1000 kg/m3.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: rhoConst :end-at: rhoConst :lineno-match:This line sets parameter :varlink:`gBaro`, the acceleration due to gravity g (in the free surface terms in :eq:`baro_model_eq_u` and :eq:`baro_model_eq_v`), to 9.81 m/s2. 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... literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: gBaro :end-at: gBaro :lineno-match:These lines set parameters :varlink:`rigidLid` and :varlink:`implicitFreeSurface` in order to suppress the rigid lid formulation of the surface pressure inverter and activate the implicit free surface formulation.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: rigidLid :end-at: implicitFreeSurface :lineno-match:This line sets parameter :varlink:`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 :filelink:`data <verification/tutorial_barotropic_gyre/input/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 :numref:`barotropic_gyre_solution`. 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... literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: momAdvection :end-at: momAdvection :lineno-match:These lines set parameters :varlink:`tempStepping` and :varlink:`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.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: tempStepping :end-at: saltStepping :lineno-match:
The first line sets the tolerance (parameter :varlink:`cg2dTargetResidual`) that the 2-D conjugate gradient solver, the iterative method used in the pressure method algorithm, will use to test for convergence. The second line sets parameter :varlink:`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.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: cg2dTargetResidual :end-at: cg2dMaxIters :lineno-match:
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 :varlink:`nIter0` is non-zero, the model would require appropriate pickup files (i.e., restart files) in order to continue integration of an existing run.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: nIter0 :end-at: nIter0 :lineno-match:This line sets parameter :varlink:`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 (:numref:`code_testing_protocols`). To integrate the solution to near steady state, uncomment the line further down where we set the value to 77760 time steps. When you make this change, be sure to also uncomment the next line that sets :varlink:`monitorFreq` (see below).
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: nTimeSteps=10 :end-at: nTimeSteps=10 :lineno-match:.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: for longer :end-at: Freq=864 :lineno-match:This line sets parameter :varlink:`deltaT`, the timestep used in stepping forward the model, to 1200 seconds. In combination with the larger value of :varlink:`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.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: deltaT :end-at: deltaT :lineno-match:
These lines control the frequency at which restart (a.k.a. pickup) files are dumped by MITgcm. Here the value of :varlink:`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 :numref:`other_output`) in order to save it as a labelled, permanent, pickup state. The value of :varlink:`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
ckptAandckptBin 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... literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: pChkptFreq :end-at: chkptFreq :lineno-match:This line sets parameter :varlink:`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.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: dumpFreq :end-at: dumpFreq :lineno-match:These lines are set to dump monitor output (see :numref:`pkg_monitor`) 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 :varlink:`monitorFreq` is set to 864,000 seconds, i.e., output every 10 days. Parameter :varlink:`monitorSelect` is set to 2 here to reduce output of non-applicable quantities for this simple example.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: monitorFreq :end-at: monitorSelect :lineno-match:
This line sets parameter :varlink:`usingCartesianGrid`, which specifies that the simulation will use a Cartesian coordinate system.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: usingCartesianGrid :end-at: usingCartesianGrid :lineno-match:These lines set the horizontal grid spacing of the model grid, as vectors :varlink:`delX` and :varlink:`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 :filelink:`SIZE.h <verification/tutorial_barotropic_gyre/code/SIZE.h>`. Grid spacing is set to 20 \times 10^{3} m (=20 km).
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: delX :end-at: delY :lineno-match:The cartesian grid default origin is (0,0) so here we set the origin with parameters :varlink:`xgOrigin` and :varlink:`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.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: xgOrigin :end-at: ygOrigin :lineno-match:This line sets parameter :varlink:`delR`, the vertical grid spacing in the z-coordinate (i.e., \Delta z), to 5000 m.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: delR :end-at: delR :lineno-match:
This line sets parameter :varlink:`bathyFile`, the name of the bathymetry file. See :numref:`baro_gyre_bathy_file` for information about the file format.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: bathyFile :end-at: bathyFile :lineno-match: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 (:varlink:`zonalWindFile`) and the y-direction (:varlink:`meridWindFile`). Note, here we have left the latter parameter blank, as there is no meridional wind stress forcing in our example.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data :start-at: zonalWindFile :end-at: meridWindFile :lineno-match:
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/data.pkg
:linenos:
:caption: verification/tutorial_barotropic_gyre/input/data.pkg
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.
.. literalinclude:: ../../../verification/tutorial_barotropic_gyre/input/eedata
:linenos:
:caption: verification/tutorial_barotropic_gyre/input/eedata
This file uses standard default values (i.e., MITgcm default is single-threaded) and does not contain customizations for this experiment.
This file is a 2-D(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”). As discussed in :numref:`sec_baro_num_config`, 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 matlab program :filelink:`verification/tutorial_barotropic_gyre/input/gendata.m` was used to generate this bathymetry file (alternatively, see python equivalent :filelink:`gendata.py <verification/tutorial_barotropic_gyre/input/gendata.py>`). By default, this file is assumed to contain 32-bit (single precision) binary numbers. See :numref:`sec_mitgcm_inp_file_format` for additional information on MITgcm input data file format specifications.
Similar to file input/bathy.bin, this file is a 2-D(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 2-D map in order
to be compatible with the standard code for loading forcing fields
in MITgcm. The matlab program :filelink:`verification/tutorial_barotropic_gyre/input/gendata.m`
was used to generate this wind stress file (alternatively, see python equivalent
:filelink:`gendata.py <verification/tutorial_barotropic_gyre/input/gendata.py>`).
To run the barotropic jet variation of this tutorial example (see :numref:`baro_jet_solution`),
you will in fact need to run one of these
programs to generate the file input/windx_siny.bin.
To configure and compile the code (following the procedure described in :numref:`building_quickstart`):
cd build ../../../tools/genmake2 -mods ../code ««-of my_platform_optionFile»» make depend make cd ..
To run the model (following the procedure in :numref:`run_the_model`):
cd run ln -s ../input/* . ln -s ../build/mitgcmuv . ./mitgcmuv > output.txt
Your run’s standard output file should be similar to :filelink:`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, :filelink:`eedata <verification/tutorial_barotropic_gyre/input/eedata>`, :filelink:`SIZE.h <verification/tutorial_barotropic_gyre/code/SIZE.h>`, :filelink:`data <verification/tutorial_barotropic_gyre/input/data>`, and :filelink:`data.pkg <verification/tutorial_barotropic_gyre/input/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 :varlink:`monitorFreq` in :filelink:`data <verification/tutorial_barotropic_gyre/input/data>`. See :numref:`pkg_monitor`).
- output from the 2-D conjugate gradient solver. More specifically, statistics from the right-hand side of the elliptic equation -- for our linear free-surface, see eq. :eq:`elliptic-backward-free-surface` -- 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 :varlink:`debugLevel` in :filelink:`data <verification/tutorial_barotropic_gyre/input/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.
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 :filelink:`utils` when binary data are read.
The following output files are generated:
Grid Data: see :numref:`spatial_discret_dyn_eq` for definitions and description of the Arakawa C-grid staggering of model variables.
XC,YC- grid cell center point locationsXG,YG- locations of grid cell verticesRC,RF- vertical cell center and cell faces positionsDXC,DYC- grid cell center point separations (:numref:`hgrid-abcd` b)DXG,DYG- separation of grid cell vertices (:numref:`hgrid-abcd` a)DRC,DRF- separation of vertical cell centers and faces, respectivelyRAC,RAS,RAW,RAZ- areas of the grid “tracer cells”, “southern cells”, “western cells” and “vorticity cells”, respectively (:numref:`hgrid-abcd`)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 2-D(x,y) data except RC, RF, DRC, DRF, which are 1-D(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 2-D 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 2-D 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
:varlink:`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 :varlink:`dumpFreq`, here set to 15,552,000 seconds, which is every 12960 time steps.
We will examine the model solutions in :numref:`barotropic_gyre_solution`.
The free-surface height is a 2-D(x,y) field.
Snapshot files exist for other prognostic model variables, in particular
filenames starting with U (:varlink:`uVel`),
V (:varlink:`uVel`), T (:varlink:`theta`), and S (:varlink:`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 3-D(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 (:varlink:`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 :varlink:`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 :varlink:`ChkptFreq`
Other Model Output Data: Model output related to reference density and hydrostatic pressure, in files Rhoref, PHrefC, PHrefF, PH, and PHL,
is discussed in depth :ref:`here <phi_hyd_discussion>` in tutorial :ref:`tutorial_baroclinic_gyre` (as these data are not terribly interesting in this single-layer setup).
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 :numref:`barotropic_nl_soln`. 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) :cite:`pedlosky:87` section 5.11.
MITgcm solution to the barotropic gyre example after t= 3.0 years of model integration. Free surface height is shown in meters.
Using matlab for example, visualizing output using the :filelink:`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 install the MITgcmutils package, see :numref:`MITgcmutils`):
from MITgcmutils import mds
import numpy as np
import matplotlib.pyplot as plt
XC = mds.rdmds('XC')
YC = mds.rdmds('YC')
Eta = mds.rdmds('Eta', 77760)
plt.contourf(XC/1000, YC/1000, Eta, np.linspace(-0.02, 0.05,8), cmap='hot_r')
plt.colorbar()
plt.show()
(for a more involved example with detailed explanations how to read in model output, perform calculations using these data, and make plots, see tutorial :ref:`Baroclinic Ocean Gyre <baroclinic_gyre_solution>`)
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 :eq:`baro_model_eq_u` and :eq:`baro_model_eq_v`.
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:
\eta(x,y) = \frac{\tau_o}{\rho_c g H} \frac{f}{\beta} \left(1 - \frac{x}{L_x}\right) \pi \sin\left(\pi \frac{y}{L_y}\right)
\left[1 - e^{-x/(2\delta_m)} \left(\cos\frac{\sqrt{3}x}{2\delta_m} + \frac{1}{\sqrt{3}} \sin\frac{\sqrt{3}x}{2\delta_m} \right) \right]
where \delta_m = (A_{h} / \beta)^{\frac{1}{3}}. :numref:`lin_anal_soln` displays the MITgcm output after switching off momentum advection vs. the analytical solution to the linearized equations. Success!
Comparison of free surface height for the near-equilibrium MITgcm solution (t= 3.0 years) with momentum advection switched off (left) and the analytical equilibrium solution to the linearized equation (right).
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 :filelink:`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 :filelink:`data <verification/tutorial_barotropic_gyre/input/data>`,
replace the line zonalWindFile='windx_cosy.bin’, with zonalWindFile='windx_siny.bin’,.
MITgcm equilibrium solution to the barotropic setup with alternate sine profile of wind stress forcing, producing a barotropic jet.
The free surface solution given this forcing is shown in :numref:`baro_jet_solution`. 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.