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 OBCSpackage (pkg/obcs) 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 compiletime and runtime options of this package, but also with the code itself. In many cases it will be necessary to adapt the obcscode (in particular S/R OBCS_CALC) to the application in question; in these cases pkg/obcs (together with the pkg/rbcs, see 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 compiletime
using the
packages.conf
file by addingobcs
to itor using
genmake2
addingenable=obcs
ordisable=obcs
switchesRequired 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 timemanagement (e.g., constant in time, or cyclic with fixed frequency) is provided through S/R OBCS_FIELDS_LOAD
More sophisticated ‘realtime’ (i.e. calendar time) management is available through S/R OBCS_PRESCRIBE_READ
The latter case requires packages pkg/cal and pkg/exf to be enabled.
Parts of the OBCS code can be enabled or disabled at compiletime 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 NOT implemented 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. Runtime parameters¶
Runtime parameters are set in files data.pkg
, data.obcs
, and
data.exf
if ‘realtime’ prescription is requested
(i.e., pkg/exf enabled). These parameter files are
read in S/Rs PACKAGES_READPARMS,
OBCS_READPARMS, and
EXF_READPARMS, respectively.
Runtime parameters may be broken into three categories:
switching on/off the package at runtime
OBCS package flags and parameters
additional timing flags in
data.exf
if selected.
8.3.1.3.1. Enabling the package¶
The OBCS package is switched on at runtime by setting
useOBCS = .TRUE.
in data.pkg
.
8.3.1.3.2. 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 

0 
Nxvector of Jindices (w.r.t. Ny) of Northern OB at each Iposition (w.r.t. Nx) 

0 
Nxvector of Jindices (w.r.t. Ny) of Southern OB at each Iposition (w.r.t. Nx) 

0 
Nyvector of Iindices (w.r.t. Nx) of Eastern OB at each Jposition (w.r.t. Ny) 

0 
Nyvector of Iindices (w.r.t. Nx) of Western OB at each Jposition (w.r.t. Ny) 

FALSE 

FALSE 

FALSE 

OBCS_balanceFacN, OBCS_balanceFacS, OBCS_balanceFacE, OBCS_balanceFacW 
1 
Factor(s) determining the details of the balancing code 
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 

2000.0 
Averaging period for phase speed (seconds) 

0.45 
Maximum allowable phase speedCFL for ABII (m/s) 

0.8 
Fixed boundary phase speed (m/s) 

FALSE 

FALSE 

Sponge layer parameters 
OBCS_PARM03 

0 
sponge layer thickness (in grid points) 

0.0 
relaxation time scale at the innermost sponge layer point of a meridional OB (s) 

0.0 
relaxation time scale at the innermost sponge layer point of a zonal OB (s) 

0.0 
relaxation time scale at the outermost sponge layer point of a meridional OB (s) 

0.0 
relaxation time scale at the outermost sponge layer point of a zonal OB (s) 

Stevens parameters 
OBCS_PARM04 

0 
Relaxation time scale for temperature/salinity (s) 

TRUE 

TRUE 
8.3.1.4. Defining open boundary positions¶
There are up to 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 OB_Jnorth(Nx) /
OB_Jsouth(Nx) and an \(N_y\)dimenisonal “column”
array OB_Ieast(Ny) / 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 Cgrid. A zero (0) element in OB_I...
/ OB_J...
means there is no corresponding OB in that column/row.
By default all elements in OB_I...
/ OB_J...
are zero. For a Northern/Southern OB, the OB Vpoint is to the South/North.
For an Eastern/Western OB, the OB Upoint is to the West/East. For example
OB_Jnorth(3)=34
means that:T(3,34)
is a an OB pointU(3,34)
is a an OB pointV(3,34)
is a an OB point
OB_Jsouth(3)=1
means that:T(3,1)
is a an OB pointU(3,1)
is a an OB pointV(3,2)
is a an OB point
OB_Ieast(10)=69
means that:T(69,10)
is a an OB pointU(69,10)
is a an OB pointV(69,10)
is a an OB point
OB_Iwest(10)=1
means that:T(1,10)
is a an OB pointU(2,10)
is a an OB pointV(1,10)
is a an OB point
For convenience, negative values for OB_Jnorth / OB_Ieast refer to points relative to the
Northern/Eastern edges of the model, e.g. OB_Jnorth(3)=1
means that the point (3,Ny)
is a northern OB
and OB_Ieast(3)=5
means that the point (3,Nx5)
is an eastern OB.
8.3.1.4.1. 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:
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,
When the boundaries are in single rows or columns as in the above example, the same can be achieved with the convenient parameters OB_singleJnorth / OB_singleJsouth / OB_singleIeast / OB_singleIwest:
OB_singleIeast = 1,
OB_singleIwest = 1,
OB_singleJnorth = 1,
OB_singleJsouth = 1,
If only the first 50 grid points of the southern boundary are boundary points:
OB_Jsouth(1:50) = 50*1,
8.3.1.4.2. A more complex example¶
Open boundaries are not restricted to single rows or columns. Each OB can be distributed in different rows and columns resulting in OBs consisting of the combination of different types of open boundaries (i.e., N, S, E and W). Figure 8.6 displays such an OB located on the leftbottom corner of a domain. Note there are five boundary points defined by southern and western boundaries. In particular, there are five southern boundary (blue lines) and two western boundaries points (red lines). For the boundary displayed in Figure 8.6 and the same dimensions as in the previous example (i.e. \(120 \times 144\) grid points), the namelist looks like this:
OB_Iwest = 1*0,1*5,142*0,
OB_Jsouth = 2*3,3*2,115*0,
For an even more complicated open boundary geometry, e.g., delimiting a concave interior domain (OB_Ieast \(\leq\) OB_Iwest), one might need to also specify the interior domain through an additional input file insideOBmaskFile for the interior mask (\(=1\) inside, \(=0\) outside).
8.3.1.5. Equations and key routines¶
8.3.1.5.1. 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).
8.3.1.5.2. OBCS_CALC:¶
Toplevel 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: OB[N/S/E/W][t/s/u/v]
; e.g. for the
salinity array at the Southern OB, the array name is
OBSs. 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 timeconstant or timevarying fields (see below).
prescribed boundary fields to compute Stevens boundary conditions.
8.3.1.5.3. ORLANSKI:¶
Orlanski radiation conditions [Orl76] examples can be found in example configurations verification/dome and verification/tutorial_plume_on_slope (as described in detail in Section 4.9).
8.3.1.5.4. OBCS_PRESCRIBE_READ:¶
When useOBCSprescribe = .TRUE.
the model tries to read
temperature, salinity, u and vvelocities from files specified in the
runtime parameters OB[N/S/E/W][t/s/u/v]File
. These files are
the usual IEEE, bigendian 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 nonlinear 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 nonhydrostatic 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 filesOB[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 external_fields_load.F or as done in pkg/exf,
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
pkg/exf is used, the time levels are
controlled for each boundary separately in the same way as the
pkg/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 OBCSWstartdate1, OBCSWstartdate2 and
OBCSWperiod. Seaice boundary values are controlled
separately with siobWstartdate1, siobWstartdate2
and siobWperiod. When pkg/exf
is not used the time levels are controlled by the runtime flags
externForcingPeriod and externForcingCycle in
data
; see verification/exp4/input/data for an example.
8.3.1.5.5. OBCS_CALC_STEVENS:¶
The boundary conditions following [Ste90] 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 equation (2.12)). (This velocity is not available at the beginning of the next time step \(n+1\), when S/Rs OBCS_CALC and 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 (either TrelaxStevens or 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 [Ste90] 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 nonlinear 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 nonzero 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.
8.3.1.5.6. OBCS_BALANCE_FLOW:¶
When turned on (CPP option ALLOW_OBCS_BALANCE defined in
OBCS_OPTIONS.h and
useOBCSbalance set to .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, OBCS_balanceFacS, OBCS_balanceFacE, and OBCS_balanceFacW. 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:
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, global balancing with OBCS_balanceFacE/W/N/S
\(\ge 0\)
is the preferred method.
Setting runtime parameter OBCSbalanceSurf to TRUE.
, the
surface mass flux contribution, say, from surface freshwater flux
EmPmR is included in the balancing scheme.
8.3.1.5.7. OBCS_APPLY_*:¶
8.3.1.5.8. OBCS_SPONGE:¶
The sponge layer code (turned on with CPP option ALLOW_OBCS_SPONGE and runtime parameter useOBCSsponge) adds a relaxation term to the righthandside 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
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 and Urelaxobcsinner set the relaxation time scales for the Eastern and Western boundaries, Vrelaxobcsbound and Vrelaxobcsinner for the Northern and Southern boundaries.
8.3.1.5.9. OB’s with nonlinear free surface¶
8.3.1.5.10. OB’s with sea ice¶
Simple Dirichlet boundary conditions for sea ice parameters can be specified in
anology to the ocean variables via filenames OB[N/S/E/W][a/h/sl/sn/u/v]File
(sea ice concentration, cell averaged sea ice thickness, salinity, cell
averaged snow thickness, ice drift components). With CPPflag
ALLOW_OBCS_SEAICE_SPONGE and runtime flags
useSeaiceSponge, seaiceSpongeThickness, and
[A/H/SL/SN]relaxobcs[inner/bound]
are available in analogy to the sponge
parameters for the ocean variables.
Neumann boundary conditions \(\frac{\partial\phi}{\partial{n}}=0\) for all sea ice variables can be applied with runtime flag SEAICEuseNeumannBC, which overrides the input files for the Dirichlet values.
Defining CPPflag OBCS_SEAICE_SMOOTH_EDGE allows to smooth the tracer seaice variables near the edges.
8.3.1.6. Flow chart¶
C !CALLING SEQUENCE:
C [...]
C  MAIN_DO_LOOP :: OpenAD case: Main timestepping loop routine
C  \ otherwise: just call FORWARD_STEP
C  
C/\  FORWARD_STEP :: Step forward a timestep ( AT LAST !!! )
C [...]
C/\   DO_OCEANIC_PHYS :: Oceanic (& seaice) physics computation
C/\    
C/\    OBCS_CALC :: Open boundary. package (see pkg/obcs).
C/\    
C [...]
C/\    SEAICE_MODEL :: pkg/seaice
C/\     SEAICE_DYNSOLVER :: pkg/seaice
C/\      OBCS_APPLY_UVICE :: apply uIce/vIce boudnary conditions
C/\     OBCS_ADJUST_UVICE :: (Only for OBCS_UVICE_OLD)
C/\     SEAICE_GROWTH
C/\     SEAICE_APPLY_SEAICE :: add OBCS for scalar variables
C [...]
C/\   THERMODYNAMICS :: theta, salt + tracer equations driver.
C/\     (synchronous timestepping case)
C [...]
C/\    TEMP_INTEGRATE :: Step forward Prognostic Eq for Temperature.
C/\    
C/\    SALT_INTEGRATE :: Step forward Prognostic Eq for Salinity.
C/\     same sequence of calls as in TEMP_INTEGRATE
C/\    
C/\    PTRACERS_INTEGRATE :: Integrate other tracer(s) (see pkg/ptracers).
C/\      same sequence of calls as in TEMP_INTEGRATE
C/\     OBCS_APPLY_PTRACER :: Open boundary package for pTracers
C/\    
C/\    OBCS_APPLY_TS :: Open boundary package (see pkg/obcs ).
C/\   
C [...]
C/\   
C/\   DYNAMICS :: Momentum equations driver.
C/\    
C [...]
C/\    OBCS_APPLY_UV :: Apply Open bndary Conditions to provisional U,V
C [...]
C/\   MOMENTUM_CORRECTION_STEP :: Finalise momentum stepping
C [...]
C/\    OBCS_APPLY_UV :: Open boundary package (see pkg/obcs).
8.3.1.7. OBCS diagnostics¶
Diagnostics output is available via the diagnostics package (see Section 9.1). Currently there are no OBCSspecific diagnostics available.
8.3.1.8. Experiments and tutorials that use obcs¶
In the directory verification the following experiments use pkg/obcs:
exp4: box with 4 open boundaries, simulating flow over a Gaussian bump based on also tests Stevensboundary conditions;
dome: based on the project “Dynamics of Overflow Mixing and Entrainment” uses OrlanskiBCs;
internal_wave: uses a heavily modified S/R OBCS_CALC
seaice_obcs: simple example who to use the seaice related code based on lab_sea;
Tutorial Gravity Plume On a Continental Slope: uses OrlanskiBCs.