8.3.1. OBCS: Open boundary conditions for regional modeling

Authors: Alistair Adcroft, Patrick Heimbach, Samar Katiwala, Martin Losch

8.3.1.1. Introduction

The OBCS-package is fundamental to regional ocean modeling with the MITgcm, but there are so many details to be considered in regional ocean modeling that this package cannot accommodate all imaginable and possible options. Therefore, for a regional simulation with very particular details it is recommended to familiarize oneself not only with the compile- and runtime-options of this package, but also with the code itself. In many cases it will be necessary to adapt the obcs-code (in particular S/R OBCS_CALC) to the application in question; in these cases the obcs-package (together with the rbcs-package Section 8.3.2) is a very useful infrastructure for implementing special regional models.

8.3.1.2. OBCS configuration and compiling

As with all MITgcm packages OBCS can be turned on or off at compile time

  • using the packages.conf file by adding obcs to it
  • or using genmake2 adding -enable=obcs or -disable=obcs switches
  • Required packages and CPP options:
    • Two alternatives are available for prescribing open boundary values, which differ in the way how OB’s are treated in time:
      • Simple time-management (e.g., constant in time, or cyclic with fixed frequency) is provided through S/R obcs_external_fields_load.
      • More sophisticated ‘real-time’ (i.e. calendar time) management is available through obcs_prescribe_read.
    • The latter case requires packages cal and exf to be enabled.

(see also Section 3.5).

Parts of the OBCS code can be enabled or disabled at compile time via CPP preprocessor flags. These options are set in OBCS_OPTIONS.h. Table 8.1 summarizes these options.

CPP option Default Description
ALLOW_OBCS_NORTH #define enable Northern OB
ALLOW_OBCS_SOUTH #define enable Southern OB
ALLOW_OBCS_EAST #define enable Eastern OB
ALLOW_OBCS_WEST #define enable Western OB
ALLOW_OBCS_PRESCRIBE #define enable code for prescribing OB’s
ALLOW_OBCS_SPONGE #undef enable sponge layer code
ALLOW_OBCS_BALANCE #define enable code for balancing transports through OB’s
ALLOW_ORLANSKI #define enable Orlanski radiation conditions at OB’s
ALLOW_OBCS_STEVENS #undef enable Stevens (1990) boundary conditions at OB’s (currently only implemented for eastern and western boundaries and NOT for ptracers)
ALLOW_OBCS_SEAICE_SPONGE #undef Include hooks to sponge layer treatment of pkg/seaice variables
ALLOW_OBCS_TIDES #undef Add tidal contributions to normal OB flow (At the moment tidal forcing is applied only to “normal” flow)

8.3.1.3. Run-time parameters

Run-time parameters are set in files data.pkg data.obcs and data.exf if ‘real-time’ prescription is requested (i.e., package exf enabled). These parameter files are read in S/R packages_readparms.F obcs_readparms.F and exf_readparms.F respectively. Run-time parameters may be broken into 3 categories:

  1. switching on/off the package at runtime
  2. OBCS package flags and parameters
  3. additional timing flags in data.exf if selected.

Enabling the package

The OBCS package is switched on at runtime by setting useOBCS = .TRUE. in data.pkg.

Package flags and parameters

Table 8.2 summarizes the runtime flags that are set in data.obcs and their default values.

Flag/parameter default Description
OB_Jnorth 0 Nx-vector of J-indices (w.r.t. Ny) of Northern OB at each I-position (w.r.t. Nx)
OB_Jsouth 0 Nx-vector of J-indices (w.r.t. Ny) of Southern OB at each I-position (w.r.t. Nx)
OB_Ieast 0 Ny-vector of I-indices (w.r.t. Nx) of Eastern OB at each J-position (w.r.t. Ny)
OB_Iwest 0 Ny-vector of I-indices (w.r.t. Nx) of Western OB at each J-position (w.r.t. Ny)
useOBCSprescribe FALSE  
useOBCSsponge FALSE  
useOBCSbalance FALSE  
OBCS_balanceFacN, OBCS_balanceFacS, OBCS_balanceFacE, OBCS_balanceFacW 1 Factor(s) determining the details of the balancing code
OBCSbalanceSurf FALSE include surface mass flux in balance
useOrlanskiNorth, useOrlanskiSouth, useOrlanskiEast, useOrlanskiWest FALSE Turn on Orlanski boundary conditions for individual boundary.
useStevensNorth, useStevensSouth, useStevensEast, useStevensWest FALSE Turn on Stevens boundary conditions for individual boundary
OBXyFile ' '

File name of OB field:

X: N(orth), S(outh), E(ast), W(est)

y: t(emperature), s(salinity), eta (sea surface height), u(-velocity), v(-velocity), w(-velocity), a (seaice area), h (sea ice thickness), sn (snow thickness), sl (sea ice salinity )
Orlanski Parameters OBCS_PARM02  
cvelTimeScale 2000.0 Averaging period for phase speed (seconds)
CMAX 0.45 Maximum allowable phase speed-CFL for AB-II (m/s)
CFIX 0.8 Fixed boundary phase speed (m/s)
useFixedCEast FALSE  
useFixedCWest FALSE  
Sponge layer parameters OBCS_PARM03  
spongeThickness 0 sponge layer thickness (in grid points)
Urelaxobcsinner 0.0 relaxation time scale at the innermost sponge layer point of a meridional OB (s)
Vrelaxobcsinner 0.0 relaxation time scale at the innermost sponge layer point of a zonal OB (s)
Urelaxobcsbound 0.0 relaxation time scale at the outermost sponge layer point of a meridional OB (s)
Vrelaxobcsbound 0.0 relaxation time scale at the outermost sponge layer point of a zonal OB (s)
Stevens parameters OBCS_PARM04  
TrelaxStevens SrelaxStevens 0 Relaxation time scale for temperature/salinity (s)
useStevensPhaseVel TRUE  
useStevensAdvection TRUE  

8.3.1.4. Defining open boundary positions

There are four open boundaries (OBs): Northern, Southern, Eastern, and Western. All OB locations are specified by their absolute meridional (Northern/Southern) or zonal (Eastern/Western) indices. Thus, for each zonal position \(i=1\ldots N_x\) a meridional index \(j\) specifies the Northern/Southern OB position, and for each meridional position \(j=1\ldots N_y\) a zonal index \(i\) specifies the Eastern/Western OB position. For Northern/Southern OB this defines an \(N_x\)-dimensional “row” array \(\tt OB\_Jnorth(Nx)\) / \(\tt OB\_Jsouth(Nx)\) and an \(N_y\)-dimenisonal “column” array \(\tt OB\_Ieast(Ny)\) / \(\tt OB\_Iwest(Ny)\). Positions determined in this way allows Northern/Southern OBs to be at variable \(j\) (or \(y\)) positions and Eastern/Western OBs at variable \(i\) (or \(x\)) positions. Here indices refer to tracer points on the C-grid. A zero (0) element in \(\tt OB\_I\ldots\) \(\tt OB\_J\ldots\) means there is no corresponding OB in that column/row. For a Northern/Southern OB, the OB V point is to the South/North. For an Eastern/Western OB, the OB U point is to the West/East. For example

OB_Jnorth(3)=34 means that:
  • T(3,34) is a an OB point
  • U(3,34) is a an OB point
  • V(3,34) is a an OB point
OB_Jsouth(3)=1 means that:
  • T(3,1) is a an OB point
  • U(3,1) is a an OB point
  • V(3,2) is a an OB point
OB_Ieast(10)=69 means that:
  • T(69,10) is a an OB point
  • U(69,10) is a an OB point
  • V(69,10) is a an OB point
OB_Iwest(10)=1 means that:
  • T(1,10) is a an OB point
  • U(2,10) is a an OB point
  • V(1,10) is a an OB point

For convenience, negative values for Jnorth/Ieast refer to points relative to the Northern/Eastern edges of the model eg. \(\tt OB\_Jnorth(3)=-1\) means that the point \(\tt (3,Ny)\) is a northern OB.

Simple examples: For a model grid with \(N_x \times N_y = 120 \times 144\) horizontal grid points with four open boundaries along the four edges of the domain the simplest way of specifying the boundary points in is:

  OB_Ieast = 144*-1,
# or OB_Ieast = 144*120,
  OB_Iwest = 144*1,
  OB_Jnorth = 120*-1,
# or OB_Jnorth = 120*144,
  OB_Jsouth = 120*1,

If only the first 50 grid points of the southern boundary are boundary points:

OB_Jsouth(1:50) = 50*1,

8.3.1.5. Equations and key routines

OBCS_READPARMS:

Set OB positions through arrays OB_Jnorth(Nx), OB_Jsouth(Nx), OB_Ieast(Ny), OB_Iwest(Ny) and runtime flags (see Table Table 8.2).

OBCS_CALC:

Top-level routine for filling values to be applied at OB for \(T,S,U,V,\eta\) into corresponding “slice” arrays \((x,z)\) \((y,z)\) for each OB: \(\tt OB[N/S/E/W][t/s/u/v]\); e.g. for salinity array at Southern OB, array name is \(\tt OBSt\). Values filled are either

  • constant vertical \(T,S\) profiles as specified in file data (tRef(Nr), sRef(Nr)) with zero velocities \(U,V\)
  • \(T,S,U,V\) values determined via Orlanski radiation conditions (see below)
  • prescribed time-constant or time-varying fields (see below).
  • use prescribed boundary fields to compute Stevens boundary conditions.

ORLANSKI:

Orlanski radiation conditions [orl:76] examples can be found in example configurations dome (http://www.rsmas.miami.edu/personal/tamay/DOME/dome.html) and plume_on_slope.

OBCS_PRESCRIBE_READ:

When useOBCSprescribe = .TRUE. the model tries to read temperature, salinity, u- and v-velocities from files specified in the runtime parameters OB[N/S/E/W][t/s/u/v]File. These files are the usual IEEE, big-endian files with dimensions of a section along an open boundary:

  • For North/South boundary files the dimensions are \((N_x\times N_r\times\mbox{time levels})\), for East/West boundary files the dimensions are \((N_y\times N_r\times\mbox{time levels})\).
  • If a non-linear free surface is used (Section 2.10.2), additional files OB[N/S/E/W]etaFile for the sea surface height $eta$ with dimension \((N_{x/y}\times\mbox{time levels})\) may be specified.
  • If non-hydrostatic dynamics are used (Section 2.9), additional files OB[N/S/E/W]wFile for the vertical velocity $w$ with dimensions \((N_{x/y}\times N_r\times\mbox{time levels})\) can be specified.
  • If useSEAICE=.TRUE. then additional files OB[N/S/E/W][a,h,sl,sn,uice,vice] for sea ice area, thickness (HEFF), seaice salinity, snow and ice velocities \((N_{x/y}\times\mbox{time levels})\) can be specified.

As in S/R external_fields_load or the exf-package, the code reads two time levels for each variable, e.g., OBNu0 and OBNu1, and interpolates linearly between these time levels to obtain the value OBNu at the current model time (step). When the exf-package is used, the time levels are controlled for each boundary separately in the same way as the exf-fields in data.exf, namelist EXF_NML_OBCS. The runtime flags follow the above naming conventions, e.g., for the western boundary the corresponding flags are OBCWstartdate1/2 and OBCWperiod. Sea-ice boundary values are controlled separately with siobWstartdate1/2 and siobWperiod. When the exf-package is not used the time levels are controlled by the runtime flags externForcingPeriod and externForcingCycle in data; see verification/exp4 for an example.

OBCS_CALC_STEVENS:

The boundary conditions following [stevens:90] require the vertically averaged normal velocity (originally specified as a stream function along the open boundary) \(\bar{u}_{ob}\) and the tracer fields \(\chi_{ob}\) (note: passive tracers are currently not implemented and the code stops when package ptracers is used together with this option). Currently the code vertically averages the normal velocity as specified in OB[E,W]u or OB[N,S]v. From these prescribed values the code computes the boundary values for the next timestep \(n+1\) as follows (as an example, we use the notation for an eastern or western boundary):

  • \(u^{n+1}(y,z) = \bar{u}_{ob}(y) + (u')^{n}(y,z)\) where \((u')^{n}\) is the deviation from the vertically averaged velocity at timestep \(n\) on the boundary. \((u')^{n}\) is computed in the previous time step \(n\) from the intermediate velocity \(u^*\) prior to the correction step (see Section 2.2 eq. (2.12)). (This velocity is not available at the beginning of the next time step \(n+1\), when S/R OBCS_CALC/OBCS_CALC_STEVENS are called, therefore it needs to be saved in S/R DYNAMICS by calling S/R OBCS\_SAVE\_UV\_N and also stored in a separate restart files pickup_stevens[N/S/E/W].${iteration}.data)

  • If \(u^{n+1}\) is directed into the model domain, the boudary value for tracer \(\chi\) is restored to the prescribed values:

    \[\chi^{n+1} = \chi^{n} + \frac{\Delta{t}}{\tau_\chi} (\chi_{ob} - \chi^{n})\]

    where \(\tau_\chi\) is the relaxation time scale T/SrelaxStevens. The new \(\chi^{n+1}\) is then subject to the advection by \(u^{n+1}\).

  • If \(u^{n+1}\) is directed out of the model domain, the tracer \(\chi^{n+1}\) on the boundary at timestep \(n+1\) is estimated from advection out of the domain with \(u^{n+1}+c\), where \(c\) is a phase velocity estimated as \(\frac{1}{2}\frac{\partial\chi}{\partial{t}}/\frac{\partial\chi}{\partial{x}}\). The numerical scheme is (as an example for an eastern boundary):

    \[\chi_{i_{b},j,k}^{n+1} = \chi_{i_{b},j,k}^{n} + \Delta{t} (u^{n+1}+c)_{i_{b},j,k}\frac{\chi_{i_{b},j,k}^{n} - \chi_{i_{b}-1,j,k}^{n}}{\Delta{x}_{i_{b}j}^{C}}\mbox{ if }u_{i_{b}jk}^{n+1}>0\]

    where \(i_{b}\) is the boundary index. For test purposes, the phase velocity contribution or the entire advection can be turned off by setting the corresponding parameters useStevensPhaseVel and useStevensAdvection to .FALSE..

See [stevens:90] for details. With this boundary condition specifying the exact net transport across the open boundary is simple, so that balancing the flow with (S/R OBCS_BALANCE_FLOW see next paragraph) is usually not necessary.

Special cases where the current implementation is not complete:

  • When you use the non-linear free surface option (parameter nonlinFreeSurf > 1), the current implementation just assumes that the gradient normal to the open boundary is zero (\(\frac{\partial\eta}{\partial{n}} = 0\)). Although this is inconsistent with geostrophic dynamics and the possibility to specify a non-zero tangent velocity together with Stevens BCs for normal velocities, it seems to work. Recommendation: Always specify zero tangential velocities with Stevens BCs.
  • There is no code for passive tracers, just a commented template in S/R obcs_calc_stevens. This means that passive tracers can be specified independently and are fluxed with the velocities that the Stevens BCs compute, but without the restoring term.
  • There are no specific Stevens BCs for sea ice, e.g., pkg/seaice. The model uses the default boundary conditions for the sea ice packages.

OBCS_BALANCE_FLOW:

When turned on (ALLOW_OBCS_BALANCE defined in OBCS_OPTIONS.h and useOBCSbalance=.true. in data.obcs/OBCS_PARM01), this routine balances the net flow across the open boundaries. By default the net flow across the boundaries is computed and all normal velocities on boundaries are adjusted to obtain zero net inflow.

This behavior can be controlled with the runtime flags OBCS_balanceFacN/S/E/W. The values of these flags determine how the net inflow is redistributed as small correction velocities between the individual sections. A value -1 balances an individual boundary, values \(>0\) determine the relative size of the correction. For example, the values

OBCS_balanceFacE = 1., OBCS_balanceFacW = -1., OBCS_balanceFacN = 2., OBCS_balanceFacS = 0.,

make the model

  • correct Western OBWu by substracting a uniform velocity to ensure zero net transport through the Western open boundary;
  • correct Eastern and Northern normal flow, with the Northern velocity correction two times larger than the Eastern correction, but not the Southern normal flow, to ensure that the total inflow through East, Northern, and Southern open boundary is balanced.

The old method of balancing the net flow for all sections individually can be recovered by setting all flags to -1. Then the normal velocities across each of the four boundaries are modified separately, so that the net volume transport across each boundary is zero. For example, for the western boundary at \(i=i_{b}\), the modified velocity is:

\[u(y,z) - \int_{\mbox{western boundary}}u dy dz \approx OBNu(j k) - \sum_{j k} OBNu(j k) h_{w}(i_{b} j k)\Delta{y_G(i_{b} j)}\Delta{z(k)}.\]

This also ensures a net total inflow of zero through all boundaries, but this combination of flags is not useful if you want to simulate, for example, a sector of the Southern Ocean with a strong ACC entering through the western and leaving through the eastern boundary, because the value of ‘’-1’’ for these flags will make sure that the strong inflow is removed. Clearly, gobal balancing with OBCS_balanceFacE/W/N/S \(\ge 0\) is the preferred method.

With runtime parameter OBCSbalanceSurf=.TRUE., the surface mass flux contribution, say, from surface freshwater flux EmPmR is included in the balancing scheme.

OBCS_APPLY_*:

OBCS_SPONGE:

The sponge layer code (turned on with ALLOW_OBCS_SPONGE and useOBCSsponge) adds a relaxation term to the right-hand-side of the momentum and tracer equations. The variables are relaxed towards the boundary values with a relaxation time scale that increases linearly with distance from the boundary

\[G_{\chi}^{\mbox{(sponge)}} = - \frac{\chi - [( L - \delta{L} ) \chi_{BC} + \delta{L}\chi]/L} {[(L-\delta{L})\tau_{b}+\delta{L}\tau_{i}]/L} = - \frac{\chi - [( 1 - l ) \chi_{BC} + l\chi]} {[(1-l)\tau_{b}+l\tau_{i}]}\]

where \(\chi\) is the model variable (U/V/T/S) in the interior, \(\chi_{BC}\) the boundary value, \(L\) the thickness of the sponge layer (runtime parameter spongeThickness in number of grid points), \(\delta{L}\in[0,L]\) (\(\frac{\delta{L}}{L}=l\in[0,1]\)) the distance from the boundary (also in grid points), and \(\tau_{b}\) (runtime parameters Urelaxobcsbound and Vrelaxobcsbound) and \(\tau_{i}\) (runtime parameters Urelaxobcsinner and Vrelaxobcsinner) the relaxation time scales on the boundary and at the interior termination of the sponge layer. The parameters Urelaxobcsbound/inner set the relaxation time scales for the Eastern and Western boundaries, Vrelaxobcsbound/inner for the Northern and Southern boundaries.

OB’s with nonlinear free surface

OB’s with sea ice

8.3.1.6. Flow chart

C     !CALLING SEQUENCE:
c ...

8.3.1.7. OBCS diagnostics

Diagnostics output is available via the diagnostics package (see Section 9). Available output fields are summarized below:

------------------------------------------------------
 <-Name->|Levs|grid|<--  Units   -->|<- Tile (max=80c)
------------------------------------------------------

8.3.1.8. Experiments and tutorials that use obcs

In the directory verification the following experiments use obcs: