Index: /branches/dev_1784_doc/DOC/NEMO_book.tex
===================================================================
--- /branches/dev_1784_doc/DOC/NEMO_book.tex	(revision 2283)
+++ /branches/dev_1784_doc/DOC/NEMO_book.tex	(revision 2284)
@@ -151,8 +151,7 @@
 \newcommand{\key} [1] {\textbf{key\_#1}\index{CPP keys!key\_#1}}	%key_cpp (key)
 \newcommand{\NEMO} {\textit{NEMO }}						%NEMO (nemo)
-% command to "commented out" portions of text ({} argument) or not ({#1} argument)
-\newcommand{\amtcomment}[1]{}   	% command to "commented out" portions of text or not (#1 in argument)
-\newcommand{\sgacomment}[1]{}   	% command to "commented out" portions of 
-\newcommand{\gmcomment}[1]{}   	% command to "commented out" portions of 
+\newcommand{\amtcomment}[1]{#1}   	% command to allow "commented out" portions of 
+\newcommand{\sgacomment}[1]{#1}   	% command to allow "commented out" portions of 
+\newcommand{\gmcomment}[1]{}   	% command to allow "commented out" portions of 
 %                                    				% text that span line breaks
 \newcommand{\alpbet} {\left(\alpha / \beta \right)} 	% alpha/beta  for slp computation
@@ -292,4 +291,8 @@
 \include{./TexFiles/Chapters/Chap_ZDF}			% Vertical diffusion
 
+\include{./TexFiles/Chapters/Chap_OBS}                  % Observation operator
+
+\include{./TexFiles/Chapters/Chap_ASM}                  % Assimilation incs
+
 \include{./TexFiles/Chapters/Chap_MISC}			% Miscellaneous topics
 
Index: /branches/dev_1784_doc/DOC/TexFiles/Chapters/Chap_ASM.tex
===================================================================
--- /branches/dev_1784_doc/DOC/TexFiles/Chapters/Chap_ASM.tex	(revision 2284)
+++ /branches/dev_1784_doc/DOC/TexFiles/Chapters/Chap_ASM.tex	(revision 2284)
@@ -0,0 +1,141 @@
+% ================================================================
+% Chapter Assimilation increments (ASM)
+% ================================================================
+\chapter{Apply assimilation increments (ASM)}
+\label{ASM}
+
+Authors: D. Lea, K. Mogensen, A. Weaver, M. Martin ...
+
+\minitoc
+
+
+\newpage
+$\ $\newline    % force a new line
+
+The ASM branch adds the functionality to apply increments to model variables, temperature, salinity, sea surface
+height, velocity and sea ice concentration. These are 
+read into the model from a file which may be produced by data assimilation. This code is controlled by the namelist
+\np{nam\_asminc}. There is a brief description of all the namelist options provided. To build the ASM code
+\np{key\_asminc} must be
+set.
+
+%===============================================================
+
+\subsection{Direct initialization}
+
+Direct initialization refers to the instantaneous correction
+of the model background state using the analysis increment.
+
+\subsection{Incremental Analysis Updates}
+
+Rather than updating the model state directly with the analysis increment,
+it may be preferable to introduce the increment gradually into the ocean
+model in order to minimize spurious adjustment processes. This technique
+is referred to as Incremental Analysis Updates (IAU; Bloom etal~1996).
+IAU is a common technique used with 3D assimilation methods such as 3D-Var
+or OI.
+
+With IAU, the model state trajectory in the assimilation window 
+($t_{0} \leq t_{i} \leq t_{N}$)
+is corrected by adding the analysis increments for 
+for temperature, salinity, horizontal velocity and SSH
+as additional tendency terms to the prognostic equations:
+\begin{eqnarray}
+{\bf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\bf x}^{b}(t_{0})] 
+\; + \; F_{i} \delta \tilde{\bf x}^{a} 
+\label{eq:wa_traj_iau}
+\end{eqnarray}
+where $F_{i}$ is a weighting function defined such that 
+$\sum_{i=1}^{N} F_{i}=1$. 
+To control the adjustment time of the model to the increment,
+the increment can be applied over an arbitrary sub-window,
+$t_{m} \leq t_{i} \leq t_{n}$, of the main assimilation window,
+where $t_{0} \leq t_{m} \leq t_{i}$ and $t_{i} \leq t_{n} \leq t_{N}$,
+Typically the increments are spread evenly over the full window.
+In addition, two different weighting functions have been implemented.
+The first function employs constant weights, 
+\begin{eqnarray}
+F^{(1)}_{i}
+=\left\{ \begin{array}{ll}
+   0   & 
+   {\rm if} \; \; \; t_{i} < t_{m} \\
+   1/M & 
+   {\rm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\
+   0  &
+   {\rm if} \; \; \; t_{i} > t_{n}
+  \end{array} \right. 
+\label{eq:F1_i}
+\end{eqnarray}
+where $M = m-n$.
+The second function employs peaked hat-like weights in order to give maximum 
+weight in the centre of the sub-window, with the weighting reduced 
+linearly to a small value at the window end-points.
+\begin{eqnarray}
+F^{(2)}_{i}
+=\left\{ \begin{array}{ll}
+   0   & 
+   {\rm if} \; \; \; t_{i} < t_{m} \\
+   \alpha \, i & 
+   {\rm if} \; \; \; t_{m} \leq t_{i} \leq t_{M/2} \\
+   \alpha \, (M - i +1) & 
+   {\rm if} \; \; \; t_{M/2} < t_{i} \leq t_{n} \\
+   0  &
+   {\rm if} \; \; \; t_{i} > t_{n}
+  \end{array} \right.
+\label{eq:F2_i}
+\end{eqnarray}
+where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even. 
+The weights described by Eq.~(\ref{eq:F2_i}) provide a 
+smoother transition of the analysis trajectory from one assimilation cycle 
+to the next than that described by Eq.~(\ref{eq:F1_i}).
+
+%==========================================================================================
+
+\section{Implementation details}
+\label{ASM_details}
+
+Here we show an example namelist and the header of an example assimilation increments file on the ORCA2 grid.
+
+\namdisplay{namasm}
+
+\subsection{Assimilation increments file}
+
+The header of an assimilation increments file produced using \np{ncdump -h} is shown below
+
+\begin{alltt}
+\tiny
+\begin{verbatim}
+netcdf assim_background_increments {
+dimensions:
+        x = 182 ;
+        y = 149 ;
+        z = 31 ;
+        t = UNLIMITED ; // (1 currently)
+variables:
+        float nav_lon(y, x) ;
+        float nav_lat(y, x) ;
+        float nav_lev(z) ;
+        double time_counter(t) ;
+        double time ;
+        double z_inc_dateb ;
+        double z_inc_datef ;
+        double bckint(t, z, y, x) ;
+        double bckins(t, z, y, x) ;
+        double bckinu(t, z, y, x) ;
+        double bckinv(t, z, y, x) ;
+        double bckineta(t, y, x) ;
+
+// global attributes:
+                :DOMAIN_number_total = 1 ;
+                :DOMAIN_number = 0 ;
+                :DOMAIN_dimensions_ids = 1, 2 ;
+                :DOMAIN_size_global = 182, 149 ;
+                :DOMAIN_size_local = 182, 149 ;
+                :DOMAIN_position_first = 1, 1 ;
+                :DOMAIN_position_last = 182, 149 ;
+                :DOMAIN_halo_size_start = 0, 0 ;
+                :DOMAIN_halo_size_end = 0, 0 ;
+                :DOMAIN_type = "BOX" ;
+}
+\end{verbatim}
+\end{alltt}
Index: /branches/dev_1784_doc/DOC/TexFiles/Chapters/Chap_OBS.tex
===================================================================
--- /branches/dev_1784_doc/DOC/TexFiles/Chapters/Chap_OBS.tex	(revision 2284)
+++ /branches/dev_1784_doc/DOC/TexFiles/Chapters/Chap_OBS.tex	(revision 2284)
@@ -0,0 +1,730 @@
+% ================================================================
+% Chapter observation operator (OBS)
+% ================================================================
+\chapter{Observation and model comparison (OBS)}
+\label{OBS}
+
+Authors: D. Lea, K. Mogensen, A. Weaver, M. Martin ...
+
+\minitoc
+
+
+\newpage
+$\ $\newline    % force a new line
+
+The OBS branch is a diagnostic branch which reads in observation files (profile temperature and
+salinity, sea surface temperature, sea level anomaly and sea ice concentration) and calculates the
+an interpolated model equivalent value at the observation location and nearest model timestep.
+This is code with was originally developed for use with NEMOVAR. 
+
+In the case of temperature data from moored buoys (TAO, TRITON, PIRATA) which in the
+ENACT/ENSEMBLES data-base are available as daily averaged quantities the observation operator
+averages the model temperature fields over one day before interpolating them to the observation
+location. For SST and altimeter observations, a 2D horizontal  interpolator is needed. For {\em in
+situ} profiles, a 1D vertical interpolator is needed in addition to the 2D interpolator.
+
+The resulting data is 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 \np{nam\_obs}. To build with the OBS code active \np{key\_diaobs} must be set.
+
+There is a brief description of all the namelist options provided. 
+
+
+% ================================================================
+% Example
+% ================================================================
+\section{Running the observation operator code example}
+\label{OBS_example}
+
+This section describes an example of running the observation operator code using
+profile data which can be freely downloaded. This shows how to adapt an existing
+run and build of NEMO to run the observation operator.
+
+First compile NEMO with \np{key\_diaobs} set.
+
+Next download some ENSEMBLES EN3 data from the website http://www.hadobs.org.
+You should choose observations which are valid for the period of your test run
+because the observation operator compares the model and observations for a
+matching date and time. 
+
+You will need to add the following to the namelist to run the observation
+operator on this data - replace \np{profiles\_01.nc} with the observation file you
+wish to use (or link in):
+
+\namdisplay{namobs_example}
+
+The option \np{ln\_t3d} and \np{ln\_s3d} switch on the temperature and salinity profile
+observation operator code. \np{ln\_ena} switch turns on the reading on ENACT/ENSEMBLES type
+profile data. The filename of the ENACT/ENSEMBLES data is set with enactfiles. Not this can be
+array of multiple files. The observation operator code needs to convert from observation
+latitude and longitude to model grid point. This is done using the grid searching part of the
+code. Setting \np{ln\_grid\_global} means that the code distributes the observations evenly
+between processors in a round-robin. Alternatively each processor will work with observations
+located within the model subdomain. The grid searching can be expensive, particularly for large
+numbers of observations, setting \np{ln\_grid\_search\_lookup} allows the use of a lookup
+table. This will need to generated the first time if it does not exist in the run directory
+however once produced it will significantly speed up future grid searches.
+
+The next step is viewing and working with the data. The NEMOVAR system contains utilities to
+plot the feedback files, convert and recombine the files which are available on request from
+the NEMOVAR team.
+
+\section{Technical details}
+\label{OBS_details}
+
+Here we show a more complete example namelist and also describe the observation files that may
+be used with the observation operator
+
+\namdisplay{namobs}
+
+This name list uses the "feedback" type observation file input format for profile, sea level
+anomaly and sea surface temperature data
+
+\subsection{Profile feedback type observation file header}
+
+\begin{alltt}
+\tiny
+\begin{verbatim}
+netcdf profiles_01 {
+dimensions:
+	N_OBS = 603 ;
+	N_LEVELS = 150 ;
+	N_VARS = 2 ;
+	N_QCF = 2 ;
+	N_ENTRIES = 1 ;
+	N_EXTRA = 1 ;
+	STRINGNAM = 8 ;
+	STRINGGRID = 1 ;
+	STRINGWMO = 8 ;
+	STRINGTYP = 4 ;
+	STRINGJULD = 14 ;
+variables:
+	char VARIABLES(N_VARS, STRINGNAM) ;
+		VARIABLES:long_name = "List of variables in feedback files" ;
+	char ENTRIES(N_ENTRIES, STRINGNAM) ;
+		ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
+	char EXTRA(N_EXTRA, STRINGNAM) ;
+		EXTRA:long_name = "List of extra variables" ;
+	char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
+		STATION_IDENTIFIER:long_name = "Station identifier" ;
+	char STATION_TYPE(N_OBS, STRINGTYP) ;
+		STATION_TYPE:long_name = "Code instrument type" ;
+	double LONGITUDE(N_OBS) ;
+		LONGITUDE:long_name = "Longitude" ;
+		LONGITUDE:units = "degrees_east" ;
+		LONGITUDE:_Fillvalue = 99999.f ;
+	double LATITUDE(N_OBS) ;
+		LATITUDE:long_name = "Latitude" ;
+		LATITUDE:units = "degrees_north" ;
+		LATITUDE:_Fillvalue = 99999.f ;
+	double DEPTH(N_OBS, N_LEVELS) ;
+		DEPTH:long_name = "Depth" ;
+		DEPTH:units = "metre" ;
+		DEPTH:_Fillvalue = 99999.f ;
+	int DEPTH_QC(N_OBS, N_LEVELS) ;
+		DEPTH_QC:long_name = "Quality on depth" ;
+		DEPTH_QC:Conventions = "q where q =[0,9]" ;
+		DEPTH_QC:_Fillvalue = 0 ;
+	int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
+		DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
+		DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+	double JULD(N_OBS) ;
+		JULD:long_name = "Julian day" ;
+		JULD:units = "days since JULD_REFERENCE" ;
+		JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
+		JULD:_Fillvalue = 99999.f ;
+	char JULD_REFERENCE(STRINGJULD) ;
+		JULD_REFERENCE:long_name = "Date of reference for julian days" ;
+		JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
+	int OBSERVATION_QC(N_OBS) ;
+		OBSERVATION_QC:long_name = "Quality on observation" ;
+		OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
+		OBSERVATION_QC:_Fillvalue = 0 ;
+	int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
+		OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
+		OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
+	int POSITION_QC(N_OBS) ;
+		POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
+		POSITION_QC:Conventions = "q where q =[0,9]" ;
+		POSITION_QC:_Fillvalue = 0 ;
+	int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
+		POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
+		POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		POSITION_QC_FLAGS:_Fillvalue = 0 ;
+	int JULD_QC(N_OBS) ;
+		JULD_QC:long_name = "Quality on date and time" ;
+		JULD_QC:Conventions = "q where q =[0,9]" ;
+		JULD_QC:_Fillvalue = 0 ;
+	int JULD_QC_FLAGS(N_OBS, N_QCF) ;
+		JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
+		JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		JULD_QC_FLAGS:_Fillvalue = 0 ;
+	int ORIGINAL_FILE_INDEX(N_OBS) ;
+		ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
+		ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
+	float POTM_OBS(N_OBS, N_LEVELS) ;
+		POTM_OBS:long_name = "Potential temperature" ;
+		POTM_OBS:units = "Degrees Celsius" ;
+		POTM_OBS:_Fillvalue = 99999.f ;
+	float POTM_Hx(N_OBS, N_LEVELS) ;
+		POTM_Hx:long_name = "Model interpolated potential temperature" ;
+		POTM_Hx:units = "Degrees Celsius" ;
+		POTM_Hx:_Fillvalue = 99999.f ;
+	int POTM_QC(N_OBS) ;
+		POTM_QC:long_name = "Quality on potential temperature" ;
+		POTM_QC:Conventions = "q where q =[0,9]" ;
+		POTM_QC:_Fillvalue = 0 ;
+	int POTM_QC_FLAGS(N_OBS, N_QCF) ;
+		POTM_QC_FLAGS:long_name = "Quality flags on potential temperature" ;
+		POTM_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		POTM_QC_FLAGS:_Fillvalue = 0 ;
+	int POTM_LEVEL_QC(N_OBS, N_LEVELS) ;
+		POTM_LEVEL_QC:long_name = "Quality for each level on potential temperature" ;
+		POTM_LEVEL_QC:Conventions = "q where q =[0,9]" ;
+		POTM_LEVEL_QC:_Fillvalue = 0 ;
+	int POTM_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
+		POTM_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on potential temperature" ;
+		POTM_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		POTM_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
+	int POTM_IOBSI(N_OBS) ;
+		POTM_IOBSI:long_name = "ORCA grid search I coordinate" ;
+	int POTM_IOBSJ(N_OBS) ;
+		POTM_IOBSJ:long_name = "ORCA grid search J coordinate" ;
+	int POTM_IOBSK(N_OBS, N_LEVELS) ;
+		POTM_IOBSK:long_name = "ORCA grid search K coordinate" ;
+	char POTM_GRID(STRINGGRID) ;
+		POTM_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
+	float PSAL_OBS(N_OBS, N_LEVELS) ;
+		PSAL_OBS:long_name = "Practical salinity" ;
+		PSAL_OBS:units = "PSU" ;
+		PSAL_OBS:_Fillvalue = 99999.f ;
+	float PSAL_Hx(N_OBS, N_LEVELS) ;
+		PSAL_Hx:long_name = "Model interpolated practical salinity" ;
+		PSAL_Hx:units = "PSU" ;
+		PSAL_Hx:_Fillvalue = 99999.f ;
+	int PSAL_QC(N_OBS) ;
+		PSAL_QC:long_name = "Quality on practical salinity" ;
+		PSAL_QC:Conventions = "q where q =[0,9]" ;
+		PSAL_QC:_Fillvalue = 0 ;
+	int PSAL_QC_FLAGS(N_OBS, N_QCF) ;
+		PSAL_QC_FLAGS:long_name = "Quality flags on practical salinity" ;
+		PSAL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		PSAL_QC_FLAGS:_Fillvalue = 0 ;
+	int PSAL_LEVEL_QC(N_OBS, N_LEVELS) ;
+		PSAL_LEVEL_QC:long_name = "Quality for each level on practical salinity" ;
+		PSAL_LEVEL_QC:Conventions = "q where q =[0,9]" ;
+		PSAL_LEVEL_QC:_Fillvalue = 0 ;
+	int PSAL_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
+		PSAL_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on practical salinity" ;
+		PSAL_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		PSAL_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
+	int PSAL_IOBSI(N_OBS) ;
+		PSAL_IOBSI:long_name = "ORCA grid search I coordinate" ;
+	int PSAL_IOBSJ(N_OBS) ;
+		PSAL_IOBSJ:long_name = "ORCA grid search J coordinate" ;
+	int PSAL_IOBSK(N_OBS, N_LEVELS) ;
+		PSAL_IOBSK:long_name = "ORCA grid search K coordinate" ;
+	char PSAL_GRID(STRINGGRID) ;
+		PSAL_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
+	float TEMP(N_OBS, N_LEVELS) ;
+		TEMP:long_name = "Insitu temperature" ;
+		TEMP:units = "Degrees Celsius" ;
+		TEMP:_Fillvalue = 99999.f ;
+
+// global attributes:
+		:title = "NEMO observation operator output" ;
+		:Convention = "NEMO unified observation operator output" ;
+}
+\end{verbatim}
+\end{alltt}
+
+\subsection{Sea level anomaly feedback type observation file header}
+
+\begin{alltt}
+\tiny
+\begin{verbatim}
+netcdf sla_01 {
+dimensions:
+	N_OBS = 41301 ;
+	N_LEVELS = 1 ;
+	N_VARS = 1 ;
+	N_QCF = 2 ;
+	N_ENTRIES = 1 ;
+	N_EXTRA = 1 ;
+	STRINGNAM = 8 ;
+	STRINGGRID = 1 ;
+	STRINGWMO = 8 ;
+	STRINGTYP = 4 ;
+	STRINGJULD = 14 ;
+variables:
+	char VARIABLES(N_VARS, STRINGNAM) ;
+		VARIABLES:long_name = "List of variables in feedback files" ;
+	char ENTRIES(N_ENTRIES, STRINGNAM) ;
+		ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
+	char EXTRA(N_EXTRA, STRINGNAM) ;
+		EXTRA:long_name = "List of extra variables" ;
+	char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
+		STATION_IDENTIFIER:long_name = "Station identifier" ;
+	char STATION_TYPE(N_OBS, STRINGTYP) ;
+		STATION_TYPE:long_name = "Code instrument type" ;
+	double LONGITUDE(N_OBS) ;
+		LONGITUDE:long_name = "Longitude" ;
+		LONGITUDE:units = "degrees_east" ;
+		LONGITUDE:_Fillvalue = 99999.f ;
+	double LATITUDE(N_OBS) ;
+		LATITUDE:long_name = "Latitude" ;
+		LATITUDE:units = "degrees_north" ;
+		LATITUDE:_Fillvalue = 99999.f ;
+	double DEPTH(N_OBS, N_LEVELS) ;
+		DEPTH:long_name = "Depth" ;
+		DEPTH:units = "metre" ;
+		DEPTH:_Fillvalue = 99999.f ;
+	int DEPTH_QC(N_OBS, N_LEVELS) ;
+		DEPTH_QC:long_name = "Quality on depth" ;
+		DEPTH_QC:Conventions = "q where q =[0,9]" ;
+		DEPTH_QC:_Fillvalue = 0 ;
+	int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
+		DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
+		DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+	double JULD(N_OBS) ;
+		JULD:long_name = "Julian day" ;
+		JULD:units = "days since JULD_REFERENCE" ;
+		JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
+		JULD:_Fillvalue = 99999.f ;
+	char JULD_REFERENCE(STRINGJULD) ;
+		JULD_REFERENCE:long_name = "Date of reference for julian days" ;
+		JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
+	int OBSERVATION_QC(N_OBS) ;
+		OBSERVATION_QC:long_name = "Quality on observation" ;
+		OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
+		OBSERVATION_QC:_Fillvalue = 0 ;
+	int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
+		OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
+		OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
+	int POSITION_QC(N_OBS) ;
+		POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
+		POSITION_QC:Conventions = "q where q =[0,9]" ;
+		POSITION_QC:_Fillvalue = 0 ;
+	int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
+		POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
+		POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		POSITION_QC_FLAGS:_Fillvalue = 0 ;
+	int JULD_QC(N_OBS) ;
+		JULD_QC:long_name = "Quality on date and time" ;
+		JULD_QC:Conventions = "q where q =[0,9]" ;
+		JULD_QC:_Fillvalue = 0 ;
+	int JULD_QC_FLAGS(N_OBS, N_QCF) ;
+		JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
+		JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		JULD_QC_FLAGS:_Fillvalue = 0 ;
+	int ORIGINAL_FILE_INDEX(N_OBS) ;
+		ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
+		ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
+	float SLA_OBS(N_OBS, N_LEVELS) ;
+		SLA_OBS:long_name = "Sea level anomaly" ;
+		SLA_OBS:units = "metre" ;
+		SLA_OBS:_Fillvalue = 99999.f ;
+	float SLA_Hx(N_OBS, N_LEVELS) ;
+		SLA_Hx:long_name = "Model interpolated sea level anomaly" ;
+		SLA_Hx:units = "metre" ;
+		SLA_Hx:_Fillvalue = 99999.f ;
+	int SLA_QC(N_OBS) ;
+		SLA_QC:long_name = "Quality on sea level anomaly" ;
+		SLA_QC:Conventions = "q where q =[0,9]" ;
+		SLA_QC:_Fillvalue = 0 ;
+	int SLA_QC_FLAGS(N_OBS, N_QCF) ;
+		SLA_QC_FLAGS:long_name = "Quality flags on sea level anomaly" ;
+		SLA_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		SLA_QC_FLAGS:_Fillvalue = 0 ;
+	int SLA_LEVEL_QC(N_OBS, N_LEVELS) ;
+		SLA_LEVEL_QC:long_name = "Quality for each level on sea level anomaly" ;
+		SLA_LEVEL_QC:Conventions = "q where q =[0,9]" ;
+		SLA_LEVEL_QC:_Fillvalue = 0 ;
+	int SLA_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
+		SLA_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on sea level anomaly" ;
+		SLA_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		SLA_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
+	int SLA_IOBSI(N_OBS) ;
+		SLA_IOBSI:long_name = "ORCA grid search I coordinate" ;
+	int SLA_IOBSJ(N_OBS) ;
+		SLA_IOBSJ:long_name = "ORCA grid search J coordinate" ;
+	int SLA_IOBSK(N_OBS, N_LEVELS) ;
+		SLA_IOBSK:long_name = "ORCA grid search K coordinate" ;
+	char SLA_GRID(STRINGGRID) ;
+		SLA_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
+	float MDT(N_OBS, N_LEVELS) ;
+		MDT:long_name = "Mean Dynamic Topography" ;
+		MDT:units = "metre" ;
+		MDT:_Fillvalue = 99999.f ;
+
+// global attributes:
+		:title = "NEMO observation operator output" ;
+		:Convention = "NEMO unified observation operator output" ;
+}
+\end{verbatim}
+\end{alltt}
+
+\subsection{Sea surface temperature feedback type observation file header}
+
+\begin{alltt}
+\tiny
+\begin{verbatim}
+netcdf sst_01 {
+dimensions:
+	N_OBS = 33099 ;
+	N_LEVELS = 1 ;
+	N_VARS = 1 ;
+	N_QCF = 2 ;
+	N_ENTRIES = 1 ;
+	STRINGNAM = 8 ;
+	STRINGGRID = 1 ;
+	STRINGWMO = 8 ;
+	STRINGTYP = 4 ;
+	STRINGJULD = 14 ;
+variables:
+	char VARIABLES(N_VARS, STRINGNAM) ;
+		VARIABLES:long_name = "List of variables in feedback files" ;
+	char ENTRIES(N_ENTRIES, STRINGNAM) ;
+		ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
+	char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
+		STATION_IDENTIFIER:long_name = "Station identifier" ;
+	char STATION_TYPE(N_OBS, STRINGTYP) ;
+		STATION_TYPE:long_name = "Code instrument type" ;
+	double LONGITUDE(N_OBS) ;
+		LONGITUDE:long_name = "Longitude" ;
+		LONGITUDE:units = "degrees_east" ;
+		LONGITUDE:_Fillvalue = 99999.f ;
+	double LATITUDE(N_OBS) ;
+		LATITUDE:long_name = "Latitude" ;
+		LATITUDE:units = "degrees_north" ;
+		LATITUDE:_Fillvalue = 99999.f ;
+	double DEPTH(N_OBS, N_LEVELS) ;
+		DEPTH:long_name = "Depth" ;
+		DEPTH:units = "metre" ;
+		DEPTH:_Fillvalue = 99999.f ;
+	int DEPTH_QC(N_OBS, N_LEVELS) ;
+		DEPTH_QC:long_name = "Quality on depth" ;
+		DEPTH_QC:Conventions = "q where q =[0,9]" ;
+		DEPTH_QC:_Fillvalue = 0 ;
+	int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
+		DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
+		DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+	double JULD(N_OBS) ;
+		JULD:long_name = "Julian day" ;
+		JULD:units = "days since JULD_REFERENCE" ;
+		JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
+		JULD:_Fillvalue = 99999.f ;
+	char JULD_REFERENCE(STRINGJULD) ;
+		JULD_REFERENCE:long_name = "Date of reference for julian days" ;
+		JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
+	int OBSERVATION_QC(N_OBS) ;
+		OBSERVATION_QC:long_name = "Quality on observation" ;
+		OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
+		OBSERVATION_QC:_Fillvalue = 0 ;
+	int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
+		OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
+		OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
+	int POSITION_QC(N_OBS) ;
+		POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
+		POSITION_QC:Conventions = "q where q =[0,9]" ;
+		POSITION_QC:_Fillvalue = 0 ;
+	int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
+		POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
+		POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		POSITION_QC_FLAGS:_Fillvalue = 0 ;
+	int JULD_QC(N_OBS) ;
+		JULD_QC:long_name = "Quality on date and time" ;
+		JULD_QC:Conventions = "q where q =[0,9]" ;
+		JULD_QC:_Fillvalue = 0 ;
+	int JULD_QC_FLAGS(N_OBS, N_QCF) ;
+		JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
+		JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		JULD_QC_FLAGS:_Fillvalue = 0 ;
+	int ORIGINAL_FILE_INDEX(N_OBS) ;
+		ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
+		ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
+	float SST_OBS(N_OBS, N_LEVELS) ;
+		SST_OBS:long_name = "Sea surface temperature" ;
+		SST_OBS:units = "Degree centigrade" ;
+		SST_OBS:_Fillvalue = 99999.f ;
+	float SST_Hx(N_OBS, N_LEVELS) ;
+		SST_Hx:long_name = "Model interpolated sea surface temperature" ;
+		SST_Hx:units = "Degree centigrade" ;
+		SST_Hx:_Fillvalue = 99999.f ;
+	int SST_QC(N_OBS) ;
+		SST_QC:long_name = "Quality on sea surface temperature" ;
+		SST_QC:Conventions = "q where q =[0,9]" ;
+		SST_QC:_Fillvalue = 0 ;
+	int SST_QC_FLAGS(N_OBS, N_QCF) ;
+		SST_QC_FLAGS:long_name = "Quality flags on sea surface temperature" ;
+		SST_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		SST_QC_FLAGS:_Fillvalue = 0 ;
+	int SST_LEVEL_QC(N_OBS, N_LEVELS) ;
+		SST_LEVEL_QC:long_name = "Quality for each level on sea surface temperature" ;
+		SST_LEVEL_QC:Conventions = "q where q =[0,9]" ;
+		SST_LEVEL_QC:_Fillvalue = 0 ;
+	int SST_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
+		SST_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on sea surface temperature" ;
+		SST_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
+		SST_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
+	int SST_IOBSI(N_OBS) ;
+		SST_IOBSI:long_name = "ORCA grid search I coordinate" ;
+	int SST_IOBSJ(N_OBS) ;
+		SST_IOBSJ:long_name = "ORCA grid search J coordinate" ;
+	int SST_IOBSK(N_OBS, N_LEVELS) ;
+		SST_IOBSK:long_name = "ORCA grid search K coordinate" ;
+	char SST_GRID(STRINGGRID) ;
+		SST_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
+
+// global attributes:
+		:title = "NEMO observation operator output" ;
+		:Convention = "NEMO unified observation operator output" ;
+}
+\end{verbatim}
+\end{alltt}
+
+\section{Theoretical details}
+
+\subsection{Horizontal interpolation methods}
+
+Consider an observation point ${\rm P}$ with 
+with longitude and latitude $({\lambda_{}}_{\rm P}, \phi_{\rm P})$ and the 
+four nearest neighbouring model grid points ${\rm A}$, ${\rm B}$, ${\rm C}$ 
+and ${\rm D}$ with longitude and latitude ($\lambda_{\rm A}$, $\phi_{\rm A}$),
+($\lambda_{\rm B}$, $\phi_{\rm B}$) etc.
+All horizontal interpolation methods implemented in NEMO
+estimate the value of a model variable $x$ at point $P$ as
+a weighted linear combination of the values of the model 
+variables at the grid points ${\rm A}$, ${\rm B}$ etc.:
+\begin{eqnarray}
+{x_{}}_{\rm P} & \hspace{-2mm} = \hspace{-2mm} & 
+\frac{1}{w} \left( {w_{}}_{\rm A} {x_{}}_{\rm A} + 
+                   {w_{}}_{\rm B} {x_{}}_{\rm B} + 
+                   {w_{}}_{\rm C} {x_{}}_{\rm C} + 
+                   {w_{}}_{\rm D} {x_{}}_{\rm D} \right)
+\end{eqnarray}
+where ${w_{}}_{\rm A}$, ${w_{}}_{\rm B}$ etc. are the respective weights for the 
+model field at points ${\rm A}$, ${\rm B}$ etc., and 
+$w = {w_{}}_{\rm A} + {w_{}}_{\rm B} + {w_{}}_{\rm C} + {w_{}}_{\rm D}$.
+
+Four different possibilities are available for computing the weights.
+
+\begin{enumerate}
+
+\item[1.] {\bf Great-Circle distance-weighted interpolation.} The weights
+  are computed as a function of the great-circle distance $s(P, \cdot)$ 
+  between $P$ and the model grid points $A$, $B$ etc. For example, 
+  the weight given to the field ${x_{}}_{\rm A}$ is specified as the 
+  product of the distances from ${\rm P}$ to the other points:
+  \begin{eqnarray}
+  {w_{}}_{\rm A} = s({\rm P}, {\rm B}) \, s({\rm P}, {\rm C}) \, s({\rm P}, {\rm D})
+  \nonumber
+  \end{eqnarray}
+  where 
+  \begin{eqnarray}
+   s\left ({\rm P}, {\rm M} \right ) 
+     & \hspace{-2mm} = \hspace{-2mm} & 
+      \cos^{-1} \! \left\{ 
+               \sin {\phi_{}}_{\rm P} \sin {\phi_{}}_{\rm M}
+             + \cos {\phi_{}}_{\rm P} \cos {\phi_{}}_{\rm M} 
+               \cos ({\lambda_{}}_{\rm M} - {\lambda_{}}_{\rm P}) 
+                   \right\}
+   \end{eqnarray}
+   and $M$ corresponds to $B$, $C$ or $D$.
+   A more stable form of the great-circle distance formula for
+   small distances ($x$ near 1) involves the arcsine function
+   (e.g., see p.~101 of Daley and Barker~(2001); NAVDAS Source Book):
+   \begin{eqnarray}
+   s\left( {\rm P}, {\rm M} \right) 
+     & \hspace{-2mm} = \hspace{-2mm} & 
+      \sin^{-1} \! \left\{ \sqrt{ 1 - x^2 } \right\}
+   \nonumber
+   \end{eqnarray}
+   where
+   \begin{eqnarray}
+    x & \hspace{-2mm} = \hspace{-2mm} & 
+      {a_{}}_{\rm M} {a_{}}_{\rm P} + {b_{}}_{\rm M} {b_{}}_{\rm P} + {c_{}}_{\rm M} {c_{}}_{\rm P}
+   \nonumber
+   \end{eqnarray}
+   and 
+   \begin{eqnarray}
+      {a_{}}_{\rm M} & \hspace{-2mm} = \hspace{-2mm} & \sin {\phi_{}}_{\rm M}, 
+      \nonumber \\
+      {a_{}}_{\rm P} & \hspace{-2mm} = \hspace{-2mm} & \sin {\phi_{}}_{\rm P}, 
+      \nonumber \\
+      {b_{}}_{\rm M} & \hspace{-2mm} = \hspace{-2mm} & \cos {\phi_{}}_{\rm M} \cos {\phi_{}}_{\rm M}, 
+      \nonumber \\
+      {b_{}}_{\rm P} & \hspace{-2mm} = \hspace{-2mm} & \cos {\phi_{}}_{\rm P} \cos {\phi_{}}_{\rm P}, 
+      \nonumber \\
+      {c_{}}_{\rm M} & \hspace{-2mm} = \hspace{-2mm} & \cos {\phi_{}}_{\rm M} \sin {\phi_{}}_{\rm M}, 
+      \nonumber \\
+      {c_{}}_{\rm P} & \hspace{-2mm} = \hspace{-2mm} & \cos {\phi_{}}_{\rm P} \sin {\phi_{}}_{\rm P}.
+      \nonumber
+   \nonumber
+  \end{eqnarray}
+
+\item[2.] {\bf Great-Circle distance-weighted interpolation with small angle
+  approximation.} Similar to the previous interpolation but with the
+  distance $s$ computed as
+  \begin{eqnarray}
+    s\left( {\rm P}, {\rm M} \right) 
+     & \hspace{-2mm} = \hspace{-2mm} & 
+      \sqrt{ \left( {\phi_{}}_{\rm M} - {\phi_{}}_{\rm P} \right)^{2} 
+      + \left( {\lambda_{}}_{\rm M} - {\lambda_{}}_{\rm P} \right)^{2}
+        \cos^{2} {\phi_{}}_{\rm M} }
+  \end{eqnarray}
+  where $M$ corresponds to $A$, $B$, $C$ or $D$.
+
+\item[3.] {\bf Bilinear interpolation for a regular spaced grid.} The
+  interpolation is split into two 1D interpolations in the longitude
+  and latitude directions, respectively.
+
+\item[4.] {\bf Bilinear remapping interpolation for a general grid.} An
+  iterative scheme that involves first mapping a quadrilateral cell
+  into a cell with coordinates (0,0), (1,0), (0,1) and (1,1). This
+  method is based on the Scrip interpolation package (Jones~2001).
+
+\end{enumerate}
+
+\subsection{Grid search}
+
+For many grids used by the NEMO model, such as the ORCA family, 
+the horizontal grid coordinates $i$ and $j$ are not simple functions 
+of latitude and longitude. Therefore, it is not always straightforward 
+to determine the grid points surrounding any given observational position.
+Before the interpolation can be performed, a search 
+algorithm is then required to determine the corner points of 
+the quadrilateral cell in which the observation is located.
+This is the most difficult and time consuming part of the 
+2D interpolation procedure. 
+A robust test for determining if an observation falls
+within a given quadrilateral cell is as follows. Let 
+${\rm P}({\lambda_{}}_{\rm P} ,{\phi_{}}_{\rm P} )$ denote the observation point,
+and let ${\rm A}({\lambda_{}}_{\rm A} ,{\phi_{}}_{\rm A} )$,
+${\rm B}({\lambda_{}}_{\rm B} ,{\phi_{}}_{\rm B} )$,
+${\rm C}({\lambda_{}}_{\rm C} ,{\phi_{}}_{\rm C} )$ 
+and 
+${\rm D}({\lambda_{}}_{\rm D} ,{\phi_{}}_{\rm D} )$ denote
+the bottom left, bottom right, top left and top right
+corner points of the cell, respectively. 
+To determine if P is inside 
+the cell, we verify that the cross-products 
+\begin{eqnarray}
+\begin{array}{lllll}
+{{\bf r}_{}}_{\rm PA} \times {{\bf r}_{}}_{\rm PC}
+& = & [({\lambda_{}}_{\rm A}\; -\; {\lambda_{}}_{\rm P} )
+      ({\phi_{}}_{\rm C}   \; -\; {\phi_{}}_{\rm P} )
+    - ({\lambda_{}}_{\rm C}\; -\; {\lambda_{}}_{\rm P} )
+      ({\phi_{}}_{\rm A}   \; -\; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
+{{\bf r}_{}}_{\rm PB} \times {{\bf r}_{}}_{\rm PA}
+& = & [({\lambda_{}}_{\rm B}\; -\; {\lambda_{}}_{\rm P} )
+      ({\phi_{}}_{\rm A}   \; -\; {\phi_{}}_{\rm P} )
+    - ({\lambda_{}}_{\rm A}\; -\; {\lambda_{}}_{\rm P} )
+      ({\phi_{}}_{\rm B}   \; -\; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
+{{\bf r}_{}}_{\rm PC} \times {{\bf r}_{}}_{\rm PD}
+& = & [({\lambda_{}}_{\rm C}\; -\; {\lambda_{}}_{\rm P} )
+      ({\phi_{}}_{\rm D}   \; -\; {\phi_{}}_{\rm P} )
+    - ({\lambda_{}}_{\rm D}\; -\; {\lambda_{}}_{\rm P} )
+      ({\phi_{}}_{\rm C}   \; -\; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
+{{\bf r}_{}}_{\rm PD} \times {{\bf r}_{}}_{\rm PB}
+& = & [({\lambda_{}}_{\rm D}\; -\; {\lambda_{}}_{\rm P} )
+      ({\phi_{}}_{\rm B}   \; -\; {\phi_{}}_{\rm P} )
+    - ({\lambda_{}}_{\rm B}\; -\; {\lambda_{}}_{\rm P} )
+      ({\phi_{}}_{\rm D}  \;  - \; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
+\end{array}
+\label{eq:cross}
+\end{eqnarray}
+point in the opposite direction to the unit normal 
+$\widehat{\bf k}$ (i.e., that the coefficients of 
+$\widehat{\bf k}$ are negative),
+where ${{\bf r}_{}}_{\rm PA}$, ${{\bf r}_{}}_{\rm PB}$, 
+etc. correspond to the vectors between points P and A, 
+P and B, etc.. The method used is
+similar to the method used in
+the Scripp interpolation package (Jones~2001).
+
+In order to speed up the grid search, there is the possibility to construct 
+a lookup table for a user specified resolution. This lookup
+table contains the lower and upper bounds on the $i$ and $j$ indices
+to be searched for on a regular grid. For each observation position,
+the closest point on the regular grid of this position is computed and
+the $i$ and $j$ ranges of this point searched to determine the precise
+four points surrounding the observation. 
+
+\subsection{Parallel aspects of horizontal interpolation}
+
+For horizontal interpolation, there is the basic problem that the
+observations are unevenly distributed on the globe. In numerical
+models, it is common to divide the model grid into subgrids (or
+domains) where each subgrid is executed on a single processing element
+with explicit message passing for exchange of information along the
+domain boundaries when running on a massively parallel processor (MPP)
+system. This approach is used by the NEMO ocean model (Madec~2008).
+
+For observations there is no natural distribution since the
+observations are not equally distributed on the globe. 
+Two options have been made available: 1) geographical distribution;
+and 2) round-robin.
+
+\subsubsection{Geographical distribution of observations among processors}
+
+\begin{figure}
+\begin{center}
+\includegraphics[width=10cm,height=12cm,angle=-90.]{./TexFiles/Figures/Fig_ASM_obsdist_local}
+\end{center}
+\caption{Example of the distribution of observations with the geographical
+distribution of observational data.} 
+\label{fig:obslocal}
+\end{figure}
+
+This is the simplest option in which the observations are distributed according 
+to the domain of the grid-point parallelization. Figure~\ref{fig:obslocal}
+shows an example of the distribution of the {\em in situ} data on processors 
+with a different colour for each observation
+on a given processor for a 4 $\times$ 2 decomposition with ORCA2. 
+The grid-point domain decomposition is clearly visible on the plot.
+
+The advantage of this approach is that all
+information needed for horizontal interpolation is available without
+any MPP communication. Of course, this is under the assumption that 
+we are only using a $2 \times 2$ grid-point stencil for the interpolation 
+(e.g., bilinear interpolation). For higher order interpolation schemes this
+is no longer valid. A disadvantage with the above scheme is that the number of
+observations on each processor can be very different. If the cost of
+the actual interpolation is expensive relative to the communication of
+data needed for interpolation, this could lead to load imbalance.
+
+\subsubsection{Round-robin distribution of observations among processors}
+
+\begin{figure}
+\begin{center}
+\includegraphics[width=10cm,height=12cm,angle=-90.]{./TexFiles/Figures/Fig_ASM_obsdist_global}
+\end{center}
+\caption{Example of the distribution of observations with the round-robin
+  distribution of observational data.}
+\label{fig:obsglobal}
+\end{figure}
+
+An alternative approach is to distribute the observations equally
+among processors and use message passing in order to retrieve 
+the stencil for interpolation. The simplest distribution of the observations 
+is to distribute them using a round-robin scheme. Figure~\ref{fig:obsglobal}
+shows the distribution of the {\em in situ} data on processors for the
+round-robin distribution of observations with a different colour for
+each observation on a given processor for a 4 $\times$ 2 decomposition 
+with ORCA2 for the same input data as in Fig.~\ref{fig:obslocal}.
+The observations are now clearly randomly distributed on the globe.
+In order to be able to perform horizontal interpolation in this case, 
+a subroutine has been developed that retrieves any grid points in the 
+global space.
+
+\subsection{Vertical interpolation operator}
+
+The 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.
Index: /branches/dev_1784_doc/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- /branches/dev_1784_doc/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2283)
+++ /branches/dev_1784_doc/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2284)
@@ -514,5 +514,5 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t] \label{Fig_SBC_dcy}  \begin{center}
-\includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_SBC_dcy.pdf}
+%\includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_SBC_dcy.pdf}
 \caption{Example of recontruction of the diurnal cycle variation of short wave flux  
 from daily mean values on an ORCA2 grid with a time sampling of 2~hours (from 1am to 11pm). 
Index: /branches/dev_1784_doc/DOC/TexFiles/Namelist/namasm
===================================================================
--- /branches/dev_1784_doc/DOC/TexFiles/Namelist/namasm	(revision 2284)
+++ /branches/dev_1784_doc/DOC/TexFiles/Namelist/namasm	(revision 2284)
@@ -0,0 +1,35 @@
+!-----------------------------------------------------------------------
+!       nam_asminc    assimilation increments namelist
+!-----------------------------------------------------------------------
+! ln_bkgwri   Logical switch for writing out background state 
+! ln_trjwri   Logical switch for writing out state trajectory
+! ln_trainc   Logical switch for applying tracer increments
+! ln_dyninc   Logical switch for applying velocity increments
+! ln_sshinc   Logical switch for applying SSH increments 
+! ln_asmdin   Logical switch for Direct Initialization (DI)
+! ln_asmiau   Logical switch for Incremental Analysis Updating (IAU)
+! nitbkg      Timestep of background in [0,nitend-nit000-1]
+! nitdin      Timestep of background for DI in [0,nitend-nit000-1]
+! nitiaustr   Timestep of start of IAU interval in [0,nitend-nit000-1]
+! nitiaufin   Timestep of end of IAU interval in [0,nitend-nit000-1]
+! niaufn      Type of IAU weighting function
+! nittrjfrq   Frequency of trajectory output for 4D-VAR
+! ln_salfix   Logical switch for ensuring that the sa > salfixmin
+! salfixmin   Minimum salinity after applying the increments
+&nam_asminc
+    ln_bkgwri = .true.
+    ln_trjwri = .false.
+    ln_trainc = .false.
+    ln_dyninc = .false.
+    ln_sshinc = .false.
+    ln_asmdin = .false.
+    ln_asmiau = .false.
+    nitbkg = 0
+    nitdin = 0
+    nitiaustr = 1
+    nitiaufin = 150
+    niaufn = 0
+    nittrjfrq = 0
+    ln_salfix = .false.
+    salfixmin = -9999
+/
Index: /branches/dev_1784_doc/DOC/TexFiles/Namelist/namobs
===================================================================
--- /branches/dev_1784_doc/DOC/TexFiles/Namelist/namobs	(revision 2284)
+++ /branches/dev_1784_doc/DOC/TexFiles/Namelist/namobs	(revision 2284)
@@ -0,0 +1,72 @@
+!-----------------------------------------------------------------------
+!       namobs    observation usage switch
+!-----------------------------------------------------------------------
+!
+!  ln_t3d                  Logical switch for T profile observations         
+!  ln_s3d                  Logical switch for S profile observations          
+!  ln_ena                  Logical switch for ENACT insitu data set           
+!  ln_cor                  Logical switch for Coriolis insitu data set       
+!  ln_profb                Logical switch for feedback insitu data set     
+!  ln_sla                  Logical switch for SLA observations               
+!  ln_sladt                Logical switch for AVISO SLA data              
+!  ln_slafb                Logical switch for feedback SLA data            
+!  ln_ssh                  Logical switch for SSH observations              
+!  ln_sst                  Logical switch for SST observations              
+!  ln_reysst               Logical switch for Reynolds observations       
+!  ln_ghrsst               Logical switch for GHRSST observations          
+!  ln_sstfb                Logical switch for feedback SST data          
+!  ln_sss                  Logical switch for SSS observations              
+!  ln_seaice               Logical switch for Sea Ice observations        
+!  ln_vel3d                Logical switch for velocity observations         
+!  ln_velavcur             Logical switch for velocity daily av. cur.    
+!  ln_velhrcur             Logical switch for velocity high freq. cur.   
+!  ln_velavadcp            Logical switch for velocity daily av. ADCP  
+!  ln_velhradcp            Logical switch for velocity high freq. ADCP
+!  ln_velfb                Logical switch for feedback velocity data       
+!  ln_grid_global          Global distribtion of observations         
+!  ln_grid_search_lookup   Logical switch for obs grid search w/lookup table  
+!  grid_search_file        Grid search lookup file header 
+!  enactfiles              ENACT input observation file names 
+!  coriofiles              Coriolis input observation file name  
+!  profbfiles              Profile feedback input observation file name 
+!  ln_profb_enatim         Enact feedback input time setting switch    
+!  slafilesact             Active SLA input observation file name
+!  slafilespas             Passive SLA input observation file name 
+!  slafbfiles              Feedback SLA input observation file name 
+!  sstfiles                GHRSST input observation file name       
+!  sstfbfiles              Feedback SST input observation file name 
+!  seaicefiles             Sea Ice input observation file name 
+!  velavcurfiles           Vel. cur. daily av. input file name  
+!  velhvcurfiles           Vel. cur. high freq. input file name  
+!  velavadcpfiles          Vel. ADCP daily av. input file name    
+!  velhvadcpfiles          Vel. ADCP high freq. input file name 
+!  velfbfiles              Vel. feedback input observation file name 
+!  dobsini                 Initial date in window YYYYMMDD.HHMMSS       
+!  dobsend                 Final date in window YYYYMMDD.HHMMSS         
+!  n1dint                  Type of vertical interpolation method        
+!  n2dint                  Type of horizontal interpolation method       
+!  ln_nea                  Rejection of observations near land switch    
+!  nmsshc                  MSSH correction scheme                         
+!  mdtcorr                 MDT  correction                               
+!  mdtcutoff               MDT cutoff for computed correction          
+!  ln_altbias              Logical switch for alt bias                
+!  ln_ignmis               Logical switch for ignoring missing files   
+!  endailyavtypes          ENACT daily average types                    
+ &namobs
+   ln_t3d = .true.
+   ln_s3d = .true.
+   ln_ena = .false.
+   ln_profb = .true.
+   ln_sla = .true.
+   ln_sladt = .false.
+   ln_slafb = .true.
+   ln_sst = .true.
+   ln_sstfb = .true.
+   profbfiles = 'profiles_01.nc'
+   slafbfiles = 'sla_01.nc'
+   sstfbfiles = 'sst_01.nc' 'sst_02.nc' 'sst_03.nc' 'sst_04.nc' 'sst_05.nc'
+   ln_altbias = .false.
+   ln_grid_global = .true.
+   ln_grid_search_lookup = .true.
+   ln_ignmis = .true.  
+/ 
Index: /branches/dev_1784_doc/DOC/TexFiles/Namelist/namobs_example
===================================================================
--- /branches/dev_1784_doc/DOC/TexFiles/Namelist/namobs_example	(revision 2284)
+++ /branches/dev_1784_doc/DOC/TexFiles/Namelist/namobs_example	(revision 2284)
@@ -0,0 +1,63 @@
+!-----------------------------------------------------------------------
+!       namobs    observation usage switch
+!-----------------------------------------------------------------------
+!
+!  ln_t3d                  Logical switch for T profile observations         
+!  ln_s3d                  Logical switch for S profile observations          
+!  ln_ena                  Logical switch for ENACT insitu data set           
+!  ln_cor                  Logical switch for Coriolis insitu data set       
+!  ln_profb                Logical switch for feedback insitu data set     
+!  ln_sla                  Logical switch for SLA observations               
+!  ln_sladt                Logical switch for AVISO SLA data              
+!  ln_slafb                Logical switch for feedback SLA data            
+!  ln_ssh                  Logical switch for SSH observations              
+!  ln_sst                  Logical switch for SST observations              
+!  ln_reysst               Logical switch for Reynolds observations       
+!  ln_ghrsst               Logical switch for GHRSST observations          
+!  ln_sstfb                Logical switch for feedback SST data          
+!  ln_sss                  Logical switch for SSS observations              
+!  ln_seaice               Logical switch for Sea Ice observations        
+!  ln_vel3d                Logical switch for velocity observations         
+!  ln_velavcur             Logical switch for velocity daily av. cur.    
+!  ln_velhrcur             Logical switch for velocity high freq. cur.   
+!  ln_velavadcp            Logical switch for velocity daily av. ADCP  
+!  ln_velhradcp            Logical switch for velocity high freq. ADCP
+!  ln_velfb                Logical switch for feedback velocity data       
+!  ln_grid_global          Global distribtion of observations         
+!  ln_grid_search_lookup   Logical switch for obs grid search w/lookup table  
+!  grid_search_file        Grid search lookup file header 
+!  enactfiles              ENACT input observation file names 
+!  coriofiles              Coriolis input observation file name  
+!  profbfiles              Profile feedback input observation file name 
+!  ln_profb_enatim         Enact feedback input time setting switch    
+!  slafilesact             Active SLA input observation file name
+!  slafilespas             Passive SLA input observation file name 
+!  slafbfiles              Feedback SLA input observation file name 
+!  sstfiles                GHRSST input observation file name       
+!  sstfbfiles              Feedback SST input observation file name 
+!  seaicefiles             Sea Ice input observation file name 
+!  velavcurfiles           Vel. cur. daily av. input file name  
+!  velhvcurfiles           Vel. cur. high freq. input file name  
+!  velavadcpfiles          Vel. ADCP daily av. input file name    
+!  velhvadcpfiles          Vel. ADCP high freq. input file name 
+!  velfbfiles              Vel. feedback input observation file name 
+!  dobsini                 Initial date in window YYYYMMDD.HHMMSS       
+!  dobsend                 Final date in window YYYYMMDD.HHMMSS         
+!  n1dint                  Type of vertical interpolation method        
+!  n2dint                  Type of horizontal interpolation method       
+!  ln_nea                  Rejection of observations near land switch    
+!  nmsshc                  MSSH correction scheme                         
+!  mdtcorr                 MDT  correction                               
+!  mdtcutoff               MDT cutoff for computed correction          
+!  ln_altbias              Logical switch for alt bias                
+!  ln_ignmis               Logical switch for ignoring missing files   
+!  endailyavtypes          ENACT daily average types                    
+ &namobs
+   ln_t3d = .true.
+   ln_s3d = .true.
+   ln_ena = .true.
+   enactfiles = 'profiles_01.nc'
+   ln_grid_global = .true.
+   ln_grid_search_lookup = .true.
+   ln_ignmis = .true.  
+/ 
