Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_ASM.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_ASM.tex	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_ASM.tex	(revision 3116)
@@ -15,6 +15,6 @@
 The ASM code adds the functionality to apply increments to the model variables: 
 temperature, salinity, sea surface height, velocity and sea ice concentration. 
-These are read into the model from a NetCDF file which may be produced by data
-assimilation.  The code can also output model background fields which are used
+These are read into the model from a NetCDF file which may be produced by separate data
+assimilation code.  The code can also output model background fields which are used
 as an input to data assimilation code. This is all controlled by the namelist
 \textit{nam\_asminc}.  There is a brief description of all the namelist options
Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 3116)
@@ -610,6 +610,5 @@
 documented or tested. 
 
-$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}=true, 
-\np{ln\_dynhpg\_hel}=true)
+$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}=true)
 \begin{equation} \label{Eq_dynhpg_sco}
 \left\{ \begin{aligned}
@@ -624,19 +623,21 @@
 \eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of 
 the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point 
-($e_{3w}$). The version \np{ln\_dynhpg\_hel}=true has been added by Aike 
-Beckmann and involves a redefinition of the relative position of $T$-points relative 
-to $w$-points. 
-
-$\bullet$ Weighted density Jacobian (WDJ) \citep{Song1998} (\np{ln\_dynhpg\_wdj}=true)
+($e_{3w}$). 
 
 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05} 
 (\np{ln\_dynhpg\_djc}=true)
 
-$\bullet$ Rotated axes scheme (rot) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_rot}=true)
-
-Note that expression \eqref{Eq_dynhpg_sco} is used when the variable volume 
+$\bullet$ Pressure Jacobian scheme (prj) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_prj}=true)
+
+Note that expression \eqref{Eq_dynhpg_sco} is commonly used when the variable volume 
 formulation is activated (\key{vvl}) because in that case, even with a flat bottom, 
 the coordinate surfaces are not horizontal but follow the free surface 
-\citep{Levier2007}. The other pressure gradient options are not yet available.
+\citep{Levier2007}. Only the pressure jacobian scheme (\np{ln\_dynhpg\_prj}=true) is available as an 
+alternative to the default \np{ln\_dynhpg\_sco}=true when \key{vvl} is active.  The pressure Jacobian scheme uses 
+a constrained cubic spline to reconstruct the density profile across the water column. This method
+maintains the monotonicity between the density nodes and is of a higher order than the linear
+interpolation method. The pressure can be calculated by analytical integration of the density profile and
+a pressure Jacobian method is used to solve the horizontal pressure gradient. This method should
+provide a more accurate calculation of the horizontal pressure gradient than the standard scheme.
 
 %--------------------------------------------------------------------------------------------------------------
@@ -1164,2 +1165,85 @@
 
 % ================================================================
+% Neptune effect 
+% ================================================================
+\section  [Neptune effect (\textit{dynnept})]
+                {Neptune effect (\mdl{dynnept})}
+\label{DYN_nept}
+
+The "Neptune effect" (thus named in \citep{HollowayOM86}) is a
+parameterisation of the potentially large effect of topographic form stress
+(caused by eddies) in driving the ocean circulation. Originally developed for
+low-resolution models, in which it was applied via a Laplacian (second-order)
+diffusion-like term in the momentum equation, it can also be applied in eddy
+permitting or resolving models, in which a more scale-selective bilaplacian
+(fourth-order) implementation is preferred. This mechanism has a
+significant effect on boundary currents (including undercurrents), and the
+upwelling of deep water near continental shelves.
+
+The theoretical basis for the method can be found in 
+\citep{HollowayJPO92}, including the explanation of why form stress is not
+necessarily a drag force, but may actually drive the flow. 
+\citep{HollowayJPO94} demonstrate the effects of the parameterisation in
+the GFDL-MOM model, at a horizontal resolution of about 1.8 degrees. 
+\citep{HollowayOM08} demonstrate the biharmonic version of the
+parameterisation in a global run of the POP model, with an average horizontal
+grid spacing of about 32km.
+
+The NEMO implementation is a simplified form of that supplied by
+Greg Holloway, the testing of which was described in \citep{HollowayJGR09}.
+The major simplification is that a time invariant Neptune velocity
+field is assumed.  This is computed only once, during start-up, and
+made available to the rest of the code via a module.  Vertical
+diffusive terms are also ignored, and the model topography itself
+is used, rather than a separate topographic dataset as in
+\citep{HollowayOM08}.  This implementation is only in the iso-level
+formulation, as is the case anyway for the bilaplacian operator.
+
+The velocity field is derived from a transport stream function given by:
+
+\begin{equation} \label{Eq_dynnept_sf}
+\psi = -fL^2H
+\end{equation}
+
+where $L$ is a latitude-dependant length scale given by:
+
+\begin{equation} \label{Eq_dynnept_ls}
+L = l_1 + (l_2 -l_1)\left ( {1 + \cos 2\phi \over 2 } \right )
+\end{equation}
+
+where $\phi$ is latitude and $l_1$ and $l_2$ are polar and equatorial length scales respectively.
+Neptune velocity components, $u^*$, $v^*$ are derived from the stremfunction as:
+
+\begin{equation} \label{Eq_dynnept_vel}
+u^* = -{1\over H} {\partial \psi \over \partial y}\ \ \  ,\ \ \ v^* = {1\over H} {\partial \psi \over \partial x}
+\end{equation}
+
+\smallskip
+%----------------------------------------------namdom----------------------------------------------------
+\namdisplay{namdyn_nept}
+%--------------------------------------------------------------------------------------------------------
+\smallskip
+
+The Neptune effect is enabled when \np{ln\_neptsimp}=true (default=false).
+\np{ln\_smooth\_neptvel} controls whether a scale-selective smoothing is applied
+to the Neptune effect flow field (default=false) (this smoothing method is as
+used by Holloway).  \np{rn\_tslse} and \np{rn\_tslsp} are the equatorial and
+polar values respectively of the length-scale parameter $L$ used in determining
+the Neptune stream function \eqref{Eq_dynnept_sf} and \eqref{Eq_dynnept_ls}.
+Values at intermediate latitudes are given by a cosine fit, mimicking the
+variation of the deformation radius with latitude.  The default values of 12km
+and 3km are those given in \citep{HollowayJPO94}, appropriate for a coarse
+resolution model. The finer resolution study of \citep{HollowayOM08} increased
+the values of L by a factor of $\sqrt 2$ to 17km and 4.2km, thus doubling the
+stream function for a given topography.
+
+The simple formulation for ($u^*$, $v^*$) can give unacceptably large velocities
+in shallow water, and \citep{HollowayOM08} add an offset to the depth in the
+denominator to control this problem. In this implementation we offer instead (at
+the suggestion of G. Madec) the option of ramping down the Neptune flow field to
+zero over a finite depth range. The switch \np{ln\_neptramp} activates this
+option (default=false), in which case velocities at depths greater than
+\np{rn\_htrmax} are unaltered, but ramp down linearly with depth to zero at a
+depth of \np{rn\_htrmin} (and shallower).
+
+% ================================================================
Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_MISC.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_MISC.tex	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_MISC.tex	(revision 3116)
@@ -253,4 +253,24 @@
 Note this implementation may be sensitive to the optimization level. 
 
+\subsection{MPP scalability}
+\label{MISC_mppsca}
+
+The default method of communicating values across the north-fold in distributed memory applications
+(\key{mpp\_mpi}) uses a \textsc{MPI\_ALLGATHER} function to exchange values from each processing
+region in the northern row with every other processing region in the northern row. This enables a
+global width array containing the top 4 rows to be collated on every northern row processor and then
+folded with a simple algorithm. Although conceptually simple, this "All to All" communication will
+hamper performance scalability for large numbers of northern row processors. From version 3.4
+onwards an alternative method is available which only performs direct "Peer to Peer" communications
+between each processor and its immediate "neighbours" across the fold line. This is achieved by
+using the default \textsc{MPI\_ALLGATHER} method during initialisation to help identify the "active"
+neighbours. Stored lists of these neighbours are then used in all subsequent north-fold exchanges to
+restrict exchanges to those between associated regions. The collated global width array for each
+region is thus only partially filled but is guaranteed to be set at all the locations actually
+required by each individual for the fold operation. This alternative method should give identical
+results to the default \textsc{ALLGATHER} method and is recommended for large values of \np{jpni}.
+The new method is activated by setting \np{ln\_nnogather} to be true ({\bf nammpp}). The
+reproducibility of results using the two methods should be confirmed for each new, non-reference
+configuration.
 
 % ================================================================
Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_OBS.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_OBS.tex	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_OBS.tex	(revision 3116)
@@ -13,39 +13,37 @@
 $\ $\newline    % force a new line
 
-The observation and model comparison code (OBS) reads in observation files
-(profile temperature and salinity, sea surface temperature, sea level anomaly,
-sea ice concentration, and velocity) and calculates  an interpolated model equivalent
-value at the observation location and nearest model timestep. The OBS code is
-called from \np{opa.F90} in order to initialise the model and to calculate the
-model equivalent values for observations on the 0th timestep. The code is then
-called again after each timestep from \np{step.F90}. The code was originally
-developed for use with NEMOVAR. 
-
-For all data types a 2D horizontal  interpolator is needed
-to interpolate the model fields to the observation location.
-For {\em in situ} profiles, a 1D vertical interpolator is needed in addition to
-provide model fields at the observation depths. Currently this only works in
-z-level model configurations, but is being developed to work with a
-generalised vertical coordinate system.
-Temperature data from moored buoys (TAO, TRITON, PIRATA) in the
-ENACT/ENSEMBLES data-base are available as daily averaged quantities. For this
-type of observation the
-observation operator will compare such observations to the model temperature
+The observation and model comparison code (OBS) reads in observation files (profile
+temperature and salinity, sea surface temperature, sea level anomaly, sea ice concentration,
+and velocity) and calculates  an interpolated model equivalent value at the observation
+location and nearest model timestep. The resulting data are saved in a ``feedback'' file (or
+files). The code was originally developed for use with the NEMOVAR data assimilation code, but
+can be used for validation or verification of model or  any other data assimilation system.
+
+The OBS code is called from \np{opa.F90} for model initialisation and to calculate the model
+equivalent values for observations on the 0th timestep. The code is then called again after
+each timestep from \np{step.F90}. To build with the OBS code active \key{diaobs} must be
+set.
+
+For all data types a 2D horizontal  interpolator is needed to interpolate the model fields to
+the observation location. For {\em in situ} profiles, a 1D vertical interpolator is needed in
+addition to provide model fields at the observation depths. Currently this only works in
+z-level model configurations, but is being developed to work with a generalised vertical
+coordinate system. Temperature data from moored buoys (TAO, TRITON, PIRATA) in the
+ENACT/ENSEMBLES data-base are available as daily averaged quantities. For this type of
+observation the observation operator will compare such observations to the model temperature
 fields averaged over one day. The relevant observation type may be specified in the namelist
-using \np{endailyavtypes}. Otherwise the model value from the nearest
-timestep to the observation time is used.
-
-The resulting data are saved in a ``feedback'' file (or files) which can be used
-for model validation and verification and also to provide information for data
-assimilation. This code is controlled by the namelist \textit{nam\_obs}. To
-build with the OBS code active \key{diaobs} must be set.
-
-Section~\ref{OBS_example} introduces a test example of the observation operator
-code including where to obtain data and how to setup the namelist.
-Section~\ref{OBS_details} introduces some more technical details of the
-different observation types used and also shows a more complete namelist.
-Finally section~\ref{OBS_theory} introduces some of the theoretical aspects of
-the observation operator including interpolation methods and running on multiple
-processors. 
+using \np{endailyavtypes}. Otherwise the model value from the nearest timestep to the
+observation time is used.
+
+The code is controlled by the namelist \textit{nam\_obs}. See the following sections for more
+details on setting up the namelist. 
+
+Section~\ref{OBS_example} introduces a test example of the observation operator code including
+where to obtain data and how to setup the namelist. Section~\ref{OBS_details} introduces some
+more technical details of the different observation types used and also shows a more complete
+namelist. Section~\ref{OBS_theory} introduces some of the theoretical aspects of the
+observation operator including interpolation methods and running on multiple processors.
+Section~\ref{OBS_obsutils} introduces some utilities to help working with the files produced
+by the OBS code.
 
 % ================================================================
@@ -69,5 +67,5 @@
 \item Add the following to the NEMO namelist to run the observation
 operator on this data. Set the \np{enactfiles} namelist parameter to the
-observation  file name (or link in to \np{profiles\_01\.nc}):
+observation  file name:
 \end{enumerate}
 
@@ -76,5 +74,5 @@
 %-------------------------------------------------------------------------------------------------------------
 
-The option \np{ln\_t3d} and \np{ln\_s3d} switch on the temperature and salinity
+The options \np{ln\_t3d} and \np{ln\_s3d} switch on the temperature and salinity
 profile observation operator code. The \np{ln\_ena} switch turns on the reading
 of ENACT/ENSEMBLES type profile data. The filename or array of filenames are
@@ -88,8 +86,8 @@
 Setting \np{ln\_grid\_global} means that the code distributes the observations
 evenly between processors. Alternatively each processor will work with
-observations located within the model subdomain.
-
-The NEMOVAR system contains utilities to plot the feedback files, convert and
-recombine the files. These are available on request from the NEMOVAR team.
+observations located within the model subdomain (see section~\ref{OBS_parallel}).
+
+A number of utilities are now provided to plot the feedback files, convert and
+recombine the files. These are explained in more detail in section~\ref{OBS_obsutils}.
 
 \section{Technical details}
@@ -713,4 +711,5 @@
 
 \subsection{Parallel aspects of horizontal interpolation}
+\label{OBS_parallel}
 
 For horizontal interpolation, there is the basic problem that the
@@ -779,7 +778,245 @@
 \subsection{Vertical interpolation operator}
 
-The vertical interpolation is achieved using either a cubic spline or
+Vertical interpolation is achieved using either a cubic spline or
 linear interpolation. For the cubic spline, the top and
 bottom boundary conditions for the second derivative of the 
 interpolating polynomial in the spline are set to zero.
 At the bottom boundary, this is done using the land-ocean mask.
+
+\newpage
+
+\section{Observation Utilities}
+\label{OBS_obsutils}
+
+Some tools for viewing and processing of observation and feedback files are provided in the
+NEMO repository for convenience. These include OBSTOOLS which are a collection of Fortran
+programs which are helpful to deal with feedback files. They do such tasks as observation file
+conversion, printing of file contents, some basic statistical analysis of feedback files. The
+other tool is an IDL program called dataplot which uses a graphical interface to visualise
+observations and feedback files. OBSTOOLS and dataplot are described in more detail below.  
+
+\subsection{Obstools}
+
+A series of Fortran utilities is provided with NEMO called OBSTOOLS. This are helpful in
+handling observation files and the feedback file output from the NEMO observation operator.
+The utilities are as follows
+
+\subsubsection{corio2fb}
+
+The program corio2fb converts profile observation files from the Coriolis format to the
+standard feedback format. The program is called in the following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+corio2fb.exe outputfile inputfile1 inputfile2 ...
+\end{verbatim}
+\end{alltt}
+
+\subsubsection{enact2fb}
+
+The program enact2fb converts profile observation files from the ENACT format to the standard
+feedback format. The program is called in the following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+enact2fb.exe outputfile inputfile1 inputfile2 ...
+\end{verbatim}
+\end{alltt}
+
+\subsubsection{fbcomb}
+
+The program fbcomb combines multiple feedback files produced by individual processors in an
+MPI run of NEMO into a single feedback file. The program is called in the following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbcomb.exe outputfile inputfile1 inputfile2 ...
+\end{verbatim}
+\end{alltt}
+
+\subsubsection{fbmatchup}
+
+The program fbmatchup will match observations from two feedback files. The program is called
+in the following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbmatchup.exe outputfile inputfile1 varname1 inputfile2 varname2 ...
+\end{verbatim}
+\end{alltt}
+
+
+\subsubsection{fbprint}
+
+The program fbprint will print the contents of a feedback file or files to standard output.
+Selected information can be output using optional arguments. The program is called in the
+following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbprint.exe [options] inputfile
+
+options:
+     -b            shorter output
+     -q            Select observations based on QC flags
+     -Q            Select observations based on QC flags
+     -B            Select observations based on QC flags
+     -u            unsorted
+     -s ID         select station ID  
+     -t TYPE       select observation type
+     -v NUM1-NUM2  select variable range to print by number 
+                      (default all)
+     -a NUM1-NUM2  select additional variable range to print by number 
+                      (default all)
+     -e NUM1-NUM2  select extra variable range to print by number 
+                      (default all)
+     -d            output date range
+     -D            print depths
+     -z            use zipped files
+\end{verbatim}
+\end{alltt}
+
+\subsubsection{fbsel}
+
+The program fbsel will select or subsample observations. The program is called in the
+following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbsel.exe <input filename> <output filename>
+\end{verbatim}
+\end{alltt}
+
+\subsubsection{fbstat}
+
+The program fbstat will output summary statistics in different global areas into a number of
+files. The program is called in the following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbstat.exe [-nmlev] <filenames>
+\end{verbatim}
+\end{alltt}
+
+\subsubsection{fbthin}
+
+The program fbthin will thin the data to 1 degree resolution. The code could easily be
+modified to thin to a different resolution. The program is called in the following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbthin.exe inputfile outputfile
+\end{verbatim}
+\end{alltt}
+
+\subsubsection{sla2fb}
+
+The program sla2fb will convert an AVISO SLA format file to feedback format. The program is
+called in the following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+sla2fb.exe [-s type] outputfile inputfile1 inputfile2 ...
+
+Option:
+     -s            Select altimeter data_source
+\end{verbatim}
+\end{alltt}
+
+\subsubsection{vel2fb}
+
+The program vel2fb will convert TAO/PIRATA/RAMA currents files to feedback format. The program
+is called in the following way:
+
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+vel2fb.exe outputfile inputfile1 inputfile2 ...
+\end{verbatim}
+\end{alltt}
+
+\subsection{building the obstools}
+
+To build the obstools use in the tools directory use ./maketools -n OBSTOOLS -m [ARCH].
+
+\subsection{Dataplot}
+
+An IDL program called dataplot is included which uses a graphical interface to visualise
+observations and feedback files. It is possible to zoom in, plot individual profiles and
+calculate some basic statistics. To plot some data run IDL and then:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+IDL> dataplot, "filename"
+\end{verbatim}
+\end{alltt}
+
+To read multiple files into dataplot, for example multiple feedback files from different
+processors or from different days, the easiest method is to use the spawn command to generate
+a list of files which can then be passed to dataplot.
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+IDL> spawn, 'ls profb*.nc', files
+IDL> dataplot, files
+\end{verbatim}
+\end{alltt}
+
+Fig~\ref{fig:obsdataplotmain} shows the main window which is launched when dataplot starts.
+This is split into three parts. At the top there is a menu bar which contains a variety of
+drop down menus. Areas - zooms into prespecified regions; plot - plots the data as a
+timeseries or a T-S diagram if appropriate; Find - allows data to be searched; Config - sets
+various configuration options.
+
+The middle part is a plot of the geographical location of the observations. This will plot the
+observation value, the model background value or observation minus background value depending
+on the option selected in the radio button at the bottom of the window. The plotting colour
+range can be changed by clicking on the colour bar. The title of the plot gives some basic
+information about the date range and depth range shown, the extreme values, and the mean and
+rms values. It is possible to zoom in using a drag-box. You may also zoom in or out using the
+mouse wheel.
+
+The bottom part of the window controls what is visible in the plot above. There are two bars
+which select the level range plotted (for profile data). The other bars below select the date
+range shown. The bottom of the figure allows the option to plot the mean, root mean square,
+standard deviation or mean square values. As mentioned above you can choose to plot the
+observation value, the model background value or observation minus background value. The next
+group of radio buttons selects the map projection. This can either be regular latitude
+longitude grid, or north or south polar stereographic. The next group of radio buttons will
+plot bad observations, switch to salinity and plot density for profile observations. The
+rightmost group of buttons will print the plot window as a postscript, save it as png, or exit
+from dataplot.
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}     \begin{center}
+%\includegraphics[width=10cm,height=12cm,angle=-90.]{./TexFiles/Figures/Fig_OBS_dataplot_main}
+\includegraphics[width=9cm,angle=-90.]{./TexFiles/Figures/Fig_OBS_dataplot_main}
+\caption{      \label{fig:obsdataplotmain}
+Main window of dataplot.}
+\end{center}     \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+If a profile point is clicked with the mouse button a plot of the observation and background
+values as a function of depth (Fig~\ref{fig:obsdataplotprofile}).
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}     \begin{center}
+%\includegraphics[width=10cm,height=12cm,angle=-90.]{./TexFiles/Figures/Fig_OBS_dataplot_prof}
+\includegraphics[width=7cm,angle=-90.]{./TexFiles/Figures/Fig_OBS_dataplot_prof}
+\caption{      \label{fig:obsdataplotprofile}
+Profile plot from dataplot produced by right clicking on a point in the main window.}
+\end{center}     \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+
+
+
Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 3116)
@@ -641,5 +641,5 @@
 \footnote{The \key{oasis4} exist. It activates portion of the code that are still under development.}. 
 It has been successfully used to interface \NEMO to most of the European atmospheric 
-GCM (ARPEGE, ECHAM, ECMWF, HadAM, LMDz), 
+GCM (ARPEGE, ECHAM, ECMWF, HadAM, HadGAM, LMDz), 
 as well as to \href{http://wrf-model.org/}{WRF} (Weather Research and Forecasting Model).
 
@@ -649,5 +649,16 @@
 When PISCES biogeochemical model (\key{top} and \key{pisces}) is also used in the coupled system, 
 the whole carbon cycle is computed by defining \key{cpl\_carbon\_cycle}. In this case, 
-CO$_2$ fluxes are exchanged between the atmosphere and the ice-ocean system.
+CO$_2$ fluxes will be exchanged between the atmosphere and the ice-ocean system (and need to be activated
+in namsbc{\_}cpl).
+
+The new namelist above allows control of various aspects of the coupling fields (particularly for
+vectors) and now allows for any coupling fields to have multiple sea ice categories (as required by LIM3
+and CICE).  When indicating a multi-category coupling field in namsbc{\_}cpl the number of categories will be
+determined by the number used in the sea ice model.  In some limited cases it may be possible to specify 
+single category coupling fields even when the sea ice model is running with multiple categories - in this
+case the user should examine the code to be sure the assumptions made are satisfactory.  In cases where
+this is definitely not possible the model should abort with an error message.  The new code has been tested using
+ECHAM with LIM2, and HadGAM3 with CICE but although it will compile with \key{lim3} additional minor code changes
+may be required to run using LIM3.
 
 
@@ -1001,10 +1012,36 @@
 ice-ocean fluxes, that are combined with the air-sea fluxes using the ice fraction of 
 each model cell to provide the surface ocean fluxes. Note that the activation of a 
-sea-ice model is is done by defining a CPP key (\key{lim2} or \key{lim3}). 
-The activation automatically ovewrite the read value of nn{\_}ice to its appropriate 
-value ($i.e.$ $2$ for LIM-2 and $3$ for LIM-3).
+sea-ice model is is done by defining a CPP key (\key{lim2}, \key{lim3} or \key{cice}). 
+The activation automatically overwrites the read value of nn{\_}ice to its appropriate 
+value ($i.e.$ $2$ for LIM-2, $3$ for LIM-3 or $4$ for CICE).
 \end{description}
 
 % {Description of Ice-ocean interface to be added here or in LIM 2 and 3 doc ?}
+
+\subsection   [Interface to CICE (\textit{sbcice\_cice})]
+			{Interface to CICE (\mdl{sbcice\_cice})}
+\label{SBC_cice}
+
+It is now possible to couple a global NEMO configuration (without AGRIF) to the CICE sea-ice
+model by using \key{cice}.  The CICE code can be obtained from 
+\href{http://oceans11.lanl.gov/trac/CICE/}{LANL} and the additional 'hadgem3' drivers will be required, 
+even with the latest code release.  Input grid files consistent with those used in NEMO will also be needed, 
+and CICE CPP keys \textbf{ORCA\_GRID}, \textbf{CICE\_IN\_NEMO} and \textbf{coupled} should be used (seek advice from UKMO 
+if necessary).  Currently the code is only designed to work when using the CORE forcing option for NEMO (with
+\textit{calc\_strair~=~true} and \textit{calc\_Tsfc~=~true} in the CICE name-list), or alternatively when NEMO 
+is coupled to the HadGAM3 atmosphere model (with \textit{calc\_strair~=~false} and \textit{calc\_Tsfc~=~false}).
+The code is intended to be used with \np{nn\_fsbc} set to 1 (although coupling ocean and ice less frequently 
+should work, it is possible the calculation of some of the ocean-ice fluxes needs to be modified slightly - the
+user should check that results are not significantly different to the standard case).
+
+There are two options for the technical coupling between NEMO and CICE.  The standard version allows
+complete flexibility for the domain decompositions in the individual models, but this is at the expense of global
+gather and scatter operations in the coupling which become very expensive on larger numbers of processors. The
+alternative option (using \key{nemocice\_decomp} for both NEMO and CICE) ensures that the domain decomposition is
+identical in both models (provided domain parameters are set appropriately, and 
+\textit{processor\_shape~=~square-ice} and \textit{distribution\_wght~=~block} in the CICE name-list) and allows
+much more efficient direct coupling on individual processors.  This solution scales much better although it is at 
+the expense of having more idle CICE processors in areas where there is no sea ice.
+
 
 % -------------------------------------------------------------------------------------------------------------
Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 3116)
@@ -1,4 +1,4 @@
 % ================================================================
-% Chapter Ñ Vertical Ocean Physics (ZDF)
+% Chapter  Vertical Ocean Physics (ZDF)
 % ================================================================
 \chapter{Vertical Ocean Physics (ZDF)}
@@ -563,5 +563,5 @@
 the clipping factor is of crucial importance for the entrainment depth predicted in 
 stably stratified situations, and that its value has to be chosen in accordance 
-with the algebraic model for the turbulent ßuxes. The clipping is only activated 
+with the algebraic model for the turbulent fluxes. The clipping is only activated 
 if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value.
 
@@ -1005,8 +1005,65 @@
 reduced as necessary to ensure stability; these changes are not reported.
 
+Limits on the bottom friction coefficient are not imposed if the user has elected to
+handle the bottom friction implicitly (see \S\ref{ZDF_bfr_imp}). The number of potential
+breaches of the explicit stability criterion are still reported for information purposes.
+
+% -------------------------------------------------------------------------------------------------------------
+%       Implicit Bottom Friction
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Implicit Bottom Friction (\np{ln\_bfrimp}$=$\textit{T})}
+\label{ZDF_bfr_imp}
+
+An optional implicit form of bottom friction has been implemented to improve
+model stability. We recommend this option for shelf sea and coastal ocean applications, especially 
+for split-explicit time splitting. This option can be invoked by setting \np{ln\_bfrimp} 
+to \textit{true} in the \textit{nambfr} namelist. This option requires \np{ln\_zdfexp} to be \textit{false} 
+in the \textit{namzdf} namelist. 
+
+This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp}, the 
+bottom boundary condition is implemented implicitly.
+
+\begin{equation} \label{Eq_dynzdf_bfr}
+\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk}
+	 = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}}
+\end{equation}
+
+where $mbk$ is the layer number of the bottom wet layer. superscript $n+1$ means the velocity used in the
+friction formula is to be calculated, so, it is implicit.
+
+If split-explicit time splitting is used, care must be taken to avoid the double counting of
+the bottom friction in the 2-D barotropic momentum equations. As NEMO only updates the barotropic 
+pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation, we need to remove
+the bottom friction induced by these two terms which has been included in the 3-D momentum trend 
+and update it with the latest value. On the other hand, the bottom friction contributed by the
+other terms (e.g. the advection term, viscosity term) has been included in the 3-D momentum equations
+and should not be added in the 2-D barotropic mode.
+
+The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the
+following:
+
+\begin{equation} \label{Eq_dynspg_ts_bfr1}
+\frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b}
+\left(\textbf{U}_{med}-\textbf{U}^{m-1}\right)
+\end{equation}
+\begin{equation} \label{Eq_dynspg_ts_bfr2}
+\frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+
+\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)-
+2\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right)
+\end{equation}
+
+where $\textbf{T}$ is the vertical integrated 3-D momentum trend. We assume the leap-frog time-stepping
+is used here. $\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step.
+ $c_{b}$ is the friction coefficient. $\eta$ is the sea surface level calculated in the barotropic loops
+while $\eta^{'}$ is the sea surface level used in the 3-D baroclinic mode. $\textbf{u}_{b}$ is the bottom
+layer horizontal velocity.
+
+
+
+
 % -------------------------------------------------------------------------------------------------------------
 %       Bottom Friction with split-explicit time splitting
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Bottom Friction with split-explicit time splitting}
+\subsection{Bottom Friction with split-explicit time splitting (\np{ln\_bfrimp}$=$\textit{F})}
 \label{ZDF_bfr_ts}
 
@@ -1017,7 +1074,7 @@
 {\key{dynspg\_flt}). Extra attention is required, however, when using 
 split-explicit time stepping (\key{dynspg\_ts}). In this case the free surface 
-equation is solved with a small time step \np{nn\_baro}*\np{rn\_rdt}, while the three 
-dimensional prognostic variables are solved with a longer time step that is a 
-multiple of \np{rn\_rdt}. The trend in the barotropic momentum due to bottom 
+equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro}, while the three 
+dimensional prognostic variables are solved with the longer time step 
+of \np{rn\_rdt} seconds. The trend in the barotropic momentum due to bottom 
 friction appropriate to this method is that given by the selected parameterisation 
 ($i.e.$ linear or non-linear bottom friction) computed with the evolving velocities 
@@ -1042,11 +1099,14 @@
 \end{enumerate}
 
-Note that the use of an implicit formulation
+Note that the use of an implicit formulation within the barotropic loop
 for the bottom friction trend means that any limiting of the bottom friction coefficient 
 in \mdl{dynbfr} does not adversely affect the solution when using split-explicit time 
 splitting. This is because the major contribution to bottom friction is likely to come from 
-the barotropic component which uses the unrestricted value of the coefficient.
-
-The implicit formulation takes the form:
+the barotropic component which uses the unrestricted value of the coefficient. However, if the
+limiting is thought to be having a major effect (a more likely prospect in coastal and shelf seas
+applications) then the fully implicit form of the bottom friction should be used (see \S\ref{ZDF_bfr_imp} ) 
+which can be selected by setting \np{ln\_bfrimp} $=$ \textit{true}.
+
+Otherwise, the implicit formulation takes the form:
 \begin{equation} \label{Eq_zdfbfr_implicitts}
  \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ]  
@@ -1115,5 +1175,5 @@
 The essential goal of the parameterization is to represent the momentum 
 exchange between the barotropic tides and the unrepresented internal waves 
-induced by the tidal ßow over rough topography in a stratified ocean. 
+induced by the tidal flow over rough topography in a stratified ocean. 
 In the current version of \NEMO, the map is built from the output of 
 the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}.
Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Introduction.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Introduction.tex	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Introduction.tex	(revision 3116)
@@ -63,4 +63,6 @@
 \citep{OASIS2006}. Two-way nesting is also available through an interface to the
 AGRIF package (Adaptative Grid Refinement in \textsc{Fortran}) \citep{Debreu_al_CG2008}.
+The interface code for coupling to an alternative sea ice model (CICE, \citet{Hunke2008}) is now 
+available although this is currently only designed for global domains, without the use of AGRIF.
 
 Other model characteristics are the lateral boundary conditions (chapter~\ref{LBC}).  
