--- manual/s_phys_pkgs/text/obcs.tex	2011/02/25 17:26:17	1.6
+++ manual/s_phys_pkgs/text/obcs.tex	2011/11/09 11:48:40	1.15
@@ -11,7 +11,17 @@
 \subsubsection{Introduction
 \label{sec:pkg:obcs:intro}}
 
-
+The OBCS-package is fundamental to regional ocean modelling with the
+MITgcm, but there are so many details to be considered in
+regional ocean modelling that this package cannot accomodate 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 \code{S/R OBCS\_CALC}) to the application
+in question; in these cases the obcs-package (together with the
+rbcs-package, section \ref{sec:pkg:rbcs}) is a very
+useful infrastructure for implementing special regional models.
 
 %----------------------------------------------------------------------
 
@@ -95,7 +105,7 @@
 \code{data.pkg}, \code{data.obcs}, and \code{data.exf} 
 if ``real-time'' prescription is requested 
 (i.e. package \code{exf} enabled).
-These parameter files are read in S/R
+vThese parameter files are read in S/R
 \code{packages\_readparms.F}, \code{obcs\_readparms.F}, and
 \code{exf\_readparms.F}, respectively.
 Run-time parameters may be broken into 3 categories:
@@ -143,6 +153,8 @@
            ~ \\
         useOBCSbalance & \code{.FALSE.} & 
            ~ \\
+           OBCS\_balanceFacN/S/E/W & 1 & factor(s) determining the details
+           of the balaning code \\
         useOrlanskiNorth/South/EastWest & \code{.FALSE.} & 
            turn on Orlanski boundary conditions for individual boundary\\
         useStevensNorth/South/EastWest & \code{.FALSE.} & 
@@ -232,21 +244,26 @@
 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.
-
-\begin{verbatim}
- For example
-     OB_Jnorth(3)=34  means that:
-          T( 3 ,34) is a an OB point
-          U(3:4,34) is a an OB point
-          V( 4 ,34) is a an OB point
- while
-     OB_Jsouth(3)=1  means that:
-          T( 3 ,1) is a an OB point
-          U(3:4,1) is a an OB point
-          V( 4 ,2) is a an OB point
-\end{verbatim}
-
-For convenience, negative values for Jnorth/Ieast refer to
+For example,
+\begin{tabbing}
+  \code{OB\_Jnorth(3)=34} \=  means that:  \= \\
+  \> \code{T(3,34)} \> is a an OB point  \\
+  \> \code{U(3,34)} \> is a an OB point \\
+  \> \code{V(3,34)} \> is a an OB point \\
+  \code{OB\_Jsouth(3)=1} \> means that: \\
+  \> \code{T(3,1)} \> is a an OB point \\
+  \> \code{U(3,1)} \> is a an OB point \\ 
+  \> \code{V(3,2)} \> is a an OB point \\
+  \code{OB\_Ieast(10)=69} \>  means that:  \>  \\
+  \> \code{T(69,10)} \> is a an OB point \\
+  \> \code{U(69,10)} \> is a an OB point \\
+  \> \code{V(69,10)} \> is a an OB point \\
+  \code{OB\_Iwest(10)=1} \>  means that:  \>  \\
+  \> \code{T(1,10)} \> is a an OB point \\
+  \> \code{U(2,10)} \> is a an OB point \\
+  \> \code{V(1,10)} \> is a an OB point
+\end{tabbing}
+For convenience, negative values for \code{Jnorth}/\code{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.
@@ -291,7 +308,6 @@
 use prescribed boundary fields to compute Stevens boundary conditions.
 \end{itemize}
 
-
 \paragraph{ORLANSKI:} ~ \\
 %
 Orlanski radiation conditions \citep{orl:76}, examples can be found in
@@ -317,66 +333,161 @@
 \item If non-hydrostatic dynamics are used
   (\ref{sec:non-hydrostatic}), additional files
   \code{OB[N/S/E/W]wFile} for the vertical velocity $w$ with
-  dimensions $(N_{x/y}\times N_r\times\mbox{time levels})$ may be
+  dimensions $(N_{x/y}\times N_r\times\mbox{time levels})$ can be
   specified.
 \item If \code{useSEAICE=.TRUE.} then additional files
   \code{OB[N/S/E/W][a,h,sl,sn,uice,vice]} for sea ice area, thickness
   (\code{HEFF}), seaice salinity, snow and ice velocities
-  $(N_{x/y}\times\mbox{time levels})$ may be specified.
+  $(N_{x/y}\times\mbox{time levels})$ can be specified.
 \end{itemize}
-When the \code{exf}-package is used, the time levels are controlled
-for each boundary separately in the same way as the \code{exf}-fields
-in \code{data.exf}, namelist \code{EXF\_NML\_OBCS}. The runtime flags
+As in \code{S/R external\_fields\_load} or the \code{exf}-package, the
+code reads two time levels for each variable, e.g.\ \code{OBNu0} and
+\code{OBNu1}, and interpolates linearly between these time levels to
+obtain the value \code{OBNu} at the current model time (step). When the
+\code{exf}-package is used, the time levels are controlled for each
+boundary separately in the same way as the \code{exf}-fields in
+\code{data.exf}, namelist \code{EXF\_NML\_OBCS}. The runtime flags
 follow the above naming conventions, e.g. for the western boundary the
 corresponding flags are \code{OBCWstartdate1/2} and
 \code{OBCWperiod}. Sea-ice boundary values are controlled separately
-with \code{siobWstartdate1/2} and \code{siobWperiod}.
-When the \code{exf}-package is not used, the time levels are
-controlled by the runtime flags \code{externForcingPeriod} and
-\code{externForcingCycle} in \code{data}, see \code{verification/exp4}
-for an example.
+with \code{siobWstartdate1/2} and \code{siobWperiod}.  When the
+\code{exf}-package is not used, the time levels are controlled by the
+runtime flags \code{externForcingPeriod} and \code{externForcingCycle}
+in \code{data}, see \code{verification/exp4} for an example.
 
 \paragraph{OBCS\_CALC\_STEVENS:} ~ \\
-(THE IMPLEMENTATION OF THESE BOUNDARY CONDITIONS IS NOT COMPLETE. SO
-FAR ONLY EASTERN AND WESTERN BOUNDARIES ARE SUPPORTED.) \\
+(THE IMPLEMENTATION OF THESE BOUNDARY CONDITIONS IS NOT
+COMPLETE. PASSIVE TRACERS, SEA ICE AND NON-LINEAR FREE SURFACE ARE NOT
+SUPPORTED PROPERLY.) \\ 
 The boundary conditions following \citet{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 \code{ptracers} is used together with this
 option). Currently, the code vertically averages the normal velocity
-as specified. 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):
+as specified in \code{OB[E,W]u} or \code{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):
 \begin{itemize}
-\item $u^{n+1}(y,z) = \bar{u}_{ob}(y) + u'(y,z)$, where $u_{n}'$ is the
-  deviation from the vertically averaged velocity one grid point
-  inward from the boundary.
+\item $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 \ref{sec:time_stepping}, e.g.,
+  eq.\,(\ref{eq:ustar-backward-free-surface})).
+  % and~(\ref{eq:vstar-backward-free-surface})). 
+  (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
+  \verb+pickup_stevens[N/S/E/W].${iteration}.data+)
+%  Define CPP-flag OBCS\_STEVENS\_USE\_INTERIOR\_VELOCITY to use the
+%  velocity one grid point inward from the boundary. 
 \item 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 \texttt{T/SrelaxStevens}.
-\item If $u^{n+1}$ is directed out of the model domain, the tracer is
-  advected 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}}$.
+  scale \texttt{T/SrelaxStevens}. The new $\chi^{n+1}$ is then subject
+  to the advection by $u^{n+1}$.
+\item 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},j,k}^{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
+  advection can be turned off by setting the corresponding parameters
   \texttt{useStevensPhaseVel} and \texttt{useStevensAdvection} to
-  \texttt{.FALSE.}.\end{itemize} See \citet{stevens:90} for details.
+  \texttt{.FALSE.}.
+\end{itemize} 
+See \citet{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.
+
+\paragraph{OBCS\_BALANCE\_FLOW:} ~ \\
+%
+When turned on (\code{ALLOW\_OBCS\_BALANCE}
+defined in \code{OBCS\_OPTIONS.h} and \code{useOBCSbalance=.true.} in
+\code{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
+\code{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 ``\code{-1}'' balances an
+individual boundary, values $>0$ determine the relative size of the
+correction. For example, the values
+\begin{tabbing}
+ \code{OBCS\_balanceFacE}\code{ = 1.,} \\
+ \code{OBCS\_balanceFacW}\code{ = -1.,} \\
+ \code{OBCS\_balanceFacN}\code{ = 2.,} \\
+ \code{OBCS\_balanceFacS}\code{ = 0.,}
+\end{tabbing}
+make the model
+\begin{itemize}
+\item correct Western \code{OBWu} by substracting a uniform velocity to
+ensure zero net transport through the Western open boundary;
+\item correct Eastern and Northern normal flow, with the Northern
+  velocity correction two times larger than the Eastern correction, but
+  \emph{not} the Southern normal flow, to ensure that the total inflow through
+  East, Northern, and Southern open boundary is balanced.
+\end{itemize}
 
-\paragraph{OBCS\_BALANCE} ~ \\
-%
-~
+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 \emph{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 \emph{not} useful if you want to
+simulate, say, a sector of the Southern Ocean with a strong ACC
+entering through the western and leaving through the eastern boundary,
+because the value of ``\code{-1}'' for these flags will make sure that
+the strong inflow is removed. Clearly, gobal balancing with
+\code{OBCS\_balanceFacE/W/N/S} $\ge0$ is the preferred method.
 
 \paragraph{OBCS\_APPLY\_*:} ~ \\
 ~
 
-\paragraph{OBCS\_SPONGE} Setting sponge layer characteristics \\
+\paragraph{OBCS\_SPONGE:} ~ \\
 %
-~
+The sponge layer code (turned on with \code{ALLOW\_OBCS\_SPONGE} and
+\code{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 \code{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 \code{Urelaxobcsbound} and
+\code{Vrelaxobcsbound}) and $\tau_{i}$ (runtime parameters
+\code{Urelaxobcsinner} and \code{Vrelaxobcsinner}) the relaxation time
+scales on the boundary and at the interior termination of the sponge
+layer. The parameters \code{Urelaxobcsbound/inner} set the relaxation
+time scales for the Eastern and Western boundaries,
+\code{Vrelaxobcsbound/inner} for the Northern and Southern boundaries.
 
 \paragraph{OB's with nonlinear free surface} ~ \\
 %
@@ -425,6 +536,22 @@
 %----------------------------------------------------------------------
 
 \subsubsection{Reference experiments}
+In the directory \code{verifcation}, the following experiments use
+\code{obcs}: 
+\begin{itemize}
+\item \code{exp4}: box with 4 open boundaries, simulating flow over a
+  Gaussian bump based on \citet{adcroft:97}, also tests
+  Stevens-boundary conditions;
+\item \code{dome}: based on the project ``Dynamics of Overflow Mixing
+  and Entrainment''
+  (\url{http://www.rsmas.miami.edu/personal/tamay/DOME/dome.html}), uses
+  Orlanski-BCs;
+\item \code{internal\_wave}: uses a heavily modified \code{S/R~OBCS\_CALC}
+\item \code{seaice\_obcs}: simple example who to use the sea-ice
+  related code, based on \code{lab\_sea};
+\item \code{tutorial\_plume\_on\_slope}: uses Orlanski-BCs, see also
+  section~\ref{sec:eg-gravityplume}.
+\end{itemize}
 
 
 
@@ -436,7 +563,7 @@
 \label{sec:pkg:obcs:experiments}
 
 \begin{itemize}
-\item{Ocean experiment in exp4 verification directory. }
+\item \code{tutorial\_plume\_on\_slope} (section~\ref{sec:eg-gravityplume})
 \end{itemize}
 
 
