--- manual/s_phys_pkgs/text/obcs.tex	2011/03/14 15:01:28	1.10
+++ manual/s_phys_pkgs/text/obcs.tex	2016/06/15 13:40:53	1.18
@@ -233,9 +233,9 @@
 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(Ny)$ / $\tt OB\_Jsouth(Ny)$,
+``row'' array $\tt OB\_Jnorth(Nx)$ / $\tt OB\_Jsouth(Nx)$,
 and an $N_y$-dimenisonal 
-``column'' array $\tt OB\_Ieast(Nx)$ / $\tt OB\_Iwest(Nx)$.
+``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.
@@ -268,6 +268,24 @@
 eg. $\tt OB\_Jnorth(3)=-1$  means that the point $\tt (3,Ny)$ 
 is a northern OB.
 
+\noindent\textbf{Simple examples:} For a model grid with $ N_{x}\times
+N_{y} = 120\times144$ horizontal grid points with four open boundaries
+along the four egdes of the domain, the simplest way of specifying the
+boundary points in \code{data.obcs} is:
+\begin{verbatim}
+  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,
+\end{verbatim}
+If only the first $50$ grid points of the southern boundary are
+boundary points: 
+\begin{verbatim}
+  OB_Jsouth(1:50) = 50*1,
+\end{verbatim}
+
 \noindent
 \textsf{Add special comments for case \#define NONLIN\_FRSURF,
 see obcs\_ini\_fixed.F}
@@ -279,7 +297,7 @@
 
 \paragraph{OBCS\_READPARMS:} ~ \\
 Set OB positions through arrays
-{\tt OB\_Jnorth(Ny), OB\_Jsouth(Ny), OB\_Ieast(Nx), OB\_Iwest(Nx)},
+{\tt OB\_Jnorth(Nx), OB\_Jsouth(Nx), OB\_Ieast(Ny), OB\_Iwest(Ny)},
 and runtime flags (see Table \ref{tab:pkg:obcs:runtime_flags}).
 
 \paragraph{OBCS\_CALC:} ~ \\
@@ -356,35 +374,60 @@
 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:} ~ \\
 %
@@ -400,21 +443,21 @@
 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, with the values
+correction. For example, the values
 \begin{tabbing}
- \code{OBCS\_balanceFac\_E}\=\code{ = 1.,} \\
- \code{OBCS\_balanceFac\_W}\>\code{ = -1.,} \\
- \code{OBCS\_balanceFac\_N}\>\code{ = 2.,} \\
- \code{OBCS\_balanceFac\_S}\>\code{ = 0.,}
+ \code{OBCS\_balanceFacE}\code{ = 1.,} \\
+ \code{OBCS\_balanceFacW}\code{ = -1.,} \\
+ \code{OBCS\_balanceFacN}\code{ = 2.,} \\
+ \code{OBCS\_balanceFacS}\code{ = 0.,}
 \end{tabbing}
-will make the model
+make the model
 \begin{itemize}
 \item correct Western \code{OBWu} by substracting a uniform velocity to
-ensure zero net transport through Western OB
+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 Eastern correction, but
-  not the Southern normal flow to ensure that the total inflow through
-  East, Northern, and Southern OB is balanced
+  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}
 
 The old method of balancing the net flow for all sections individually
@@ -427,13 +470,13 @@
 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 to
-make it a useful flag for preventing infinite sea-level change within
-the domain, 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.
+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\_*:} ~ \\
 ~
@@ -455,7 +498,7 @@
 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]$ ($l\in[0,1]$) the distance from the boundary (also in grid points), and
+$\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