Index: trunk/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- trunk/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2531)
+++ trunk/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2541)
@@ -63,5 +63,4 @@
 \label{SBC_general}
 
-
 The surface ocean stress is the stress exerted by the wind and the sea-ice 
 on the ocean. The two components of stress are assumed to be interpolated 
@@ -176,4 +175,260 @@
 
 % ================================================================
+%       Input Data 
+% ================================================================
+\section{Input Data generic interface}
+\label{SBC_input}
+
+A generic interface has been introduced to manage the way input data (2D or 3D fields, 
+like surface forcing or ocean T and S) are specify in \NEMO. This task is archieved by fldread.F90. 
+The module was design with four main objectives in mind: 
+\begin{enumerate}  
+\item optionally provide a time interpolation of the input data at model time-step, 
+whatever their input frequency is, and according to the different calendars available in the model.
+\item optionally provide an on-the-fly space interpolation from the native input data grid to the model grid.
+\item make the run duration independent from the period cover by the input files.
+\item provide a simple user interface and a rather simple developer interface by limiting the 
+ number of prerequisite information. 
+\end{enumerate}  
+
+As a results the user have only to fill in for each variable a structure in the namelist file 
+to defined the input data file and variable names, the frequency of the data (in hours or months), 
+whether its is climatological data or not, the period covered by the input file (one year, month, week or day), 
+and two additional parameters for on-the-fly interpolation. When adding a new input variable, 
+the developer has to add the associated structure in the namelist, read this information 
+by mirroring the namelist read in \rou{sbc\_blk\_init} for example, and simply call \rou{fld\_read} 
+to obtain the desired input field at the model time-step and grid points.
+
+The only constraints are that the input file is a NetCDF file, the file name follows a nomenclature 
+(see \S\ref{SBC_fldread}), the period it cover is one year, month, week or day, and, if on-the-fly 
+interpolation is used, a file of weights must be supplied (see \S\ref{SBC_iof}).
+
+Note that when an input data is archived on a disc which is accessible directly 
+from the workspace where the code is executed, then the use can set the \np{cn\_dir} 
+to the pathway leading to the data. By default, the data are assumed to have been 
+copied so that cn\_dir='./'.
+
+% -------------------------------------------------------------------------------------------------------------
+% Input Data specification (\mdl{fldread})
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Input Data specification (\mdl{fldread})}
+\label{SBC_fldread}
+
+The structure associated with an input variable contains the following information:
+\begin{alltt}  {{\tiny    
+\begin{verbatim}
+!  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! 
+!             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+\end{verbatim}
+}}\end{alltt} 
+where 
+\begin{description}  
+\item[File name]: the stem name of the NetCDF file to be open. 
+This stem will be completed automatically by the model, with the addition of a '.nc' at its end 
+and by date information and possibly a prefix (when using AGRIF). 
+Tab.\ref{Tab_fldread} provides the resulting file name in all possible cases according to whether 
+it is a climatological file or not, and to the open/close frequency (see below for definition). 
+
+%--------------------------------------------------TABLE--------------------------------------------------
+\begin{table}[htbp] 
+\begin{center}
+\begin{tabular}{|l|c|c|c|}
+\hline
+                         & daily or weekLLL	        & monthly                   &   yearly          \\   \hline
+clim = false	& fn\_yYYYYmMMdDD  &   fn\_yYYYYmMM   &   fn\_yYYYY  \\   \hline
+clim = true	 	   & not possible 	              &  fn\_m??.nc             &   fn                \\   \hline
+\end{tabular}
+\end{center}
+\caption{ \label{Tab_fldread}   naming nomenclature for climatological or interannual input file, 
+as a function of the Open/close frequency. The stem name is assumed to be 'fn'. 
+For weekly files, the 'LLL' corresponds to the first three letters of the first day of the week ($i.e.$ 'sun','sat','fri','thu','wed','tue','mon'). The 'YYYY', 'MM' and 'DD' should be replaced by the 
+actual year/month/day, always coded with 4 or 2 digits. Note that (1) in mpp, if the file is split 
+over each subdomain, the suffix '.nc' is replaced by '\_PPPP.nc', where 'PPPP' is the 
+process number coded with 4 digits; (2) when using AGRIF, the prefix ÔN\_Õ is added to files, 
+where 'N'  is the child grid number.}
+\end{table}
+%--------------------------------------------------------------------------------------------------------------
+  
+
+\item[Record frequency]: the frequency of the records contained in the input file. 
+Its unit is in hours if it is positive (for example 24 for daily forcing) or in months if negative 
+(for example -1 for monthly forcing or -12 for annual forcing). 
+Note that this frequency must really be an integer and not a real. 
+On some computers, seting it to '24.' can be interpreted as 240!
+
+\item[Variable name]: the name of the variable to be read in the input NetCDF file.
+
+\item[Time interpolation]: a logical to activate, or not, the time interpolation. If set to 'false', 
+the forcing will have a steplike shape remaining constant during each forcing period. 
+For example, when using a daily forcing without time interpolation, the forcing remaining 
+constant from 00h00'00'' to 23h59'59". If set to 'true', the forcing will have a broken line shape. 
+Records are assumed to be dated the middle of the forcing period. 
+For example, when using a daily forcing with time interpolation, linear interpolation will 
+be performed between mid-day of two consecutive days. 
+
+\item[Climatological forcing]: a logical to specify if a input file contains climatological forcing 
+which can be cycle in time, or an interannual forcing which will requires additional files 
+if the period covered by the simulation exceed the one of the file. See the above the file 
+naming strategy which impacts the expected name of the file to be opened. 
+
+\item[Open/close frequency]: the frequency at which forcing files must be opened/closed. 
+Four cases are coded: 'daily', 'weekLLL' (with 'LLL' the first 3 letters of the first day of the week), 
+'monthly' and 'yearly' which means the forcing files will contain data for one day, one week, 
+one month or one year. Files are assumed to contain data from the beginning of the open/close period. 
+For example, the first record of a yearly file containing daily data is Jan 1st even if the experiment 
+is not starting at the beginning of the year. 
+
+\item[Others]: 'weights filename' and 'pairing rotation' are associted with on-the-fly interpolation 
+which is described in \S\ref{SBC_iof}.
+
+\end{description}
+
+Additional remarks:\\
+(1) The time interpolation is a simple linear interpolation between two consecutive records of 
+the input data. The only tricky point is therefore to specify the date at which we need to do 
+the interpolation and the date of the records read in the input files. 
+Following \citet{Leclair_Madec_OM09}, the date of a time step is set at the middle of the 
+time step. For example, for an experiment starting at 0h00'00" with a one hour time-step, 
+a time interpolation will be performed at the following time: 0h30'00", 1h30'00", 2h30'00", etc.
+However, for forcing data related to the surface module, values are not needed at every 
+time-step but at every \np{nn\_fsbc} time-step. For example with \np{nn\_fsbc}~=~3, 
+the surface module will be called at time-steps 1, 4, 7, etc. The date used for the time interpolation 
+is thus redefined to be at the middle of \np{nn\_fsbc} time-step period. In the previous example, 
+this leads to: 1h30'00", 4h30'00", 7h30'00", etc. \\ 
+(2) For code readablility and maintenance issues, we don't take into account the NetCDF input file 
+calendar. The calendar associated with the forcing field is build according to the information 
+provided by user in the record frequency, the open/close frequency and the type of temporal interpolation. 
+For example, the first record of a yearly file containing daily data that will be interpolated in time 
+is assumed to be start Jan 1st at 12h00'00" and end Dec 31st at 12h00'00". \\
+(3) If a time interpolation is requested, the code will pick up the needed data in the previous (next) file 
+when interpolating data with the first (last) record of the open/close period. 
+For example, if the input file specifications are ''yearly, containing daily data to be interpolated in time'', 
+the values given by the code between 00h00'00" and 11h59'59" on Jan 1st will be interpolated values 
+between Dec 31st 12h00'00" and Jan 1st 12h00'00". If the forcing is climatological, Dec and Jan will 
+be keep-up from the same year. However, if the forcing is not climatological, at the end of the 
+open/close period the code will automatically close the current file and open the next one. 
+Note that, if the experiment is starting (ending) at the beginning (end) of an open/close period 
+we do accept that the previous (next) file is not existing. In this case, the time interpolation 
+will be performed between two identical values. For example, when starting an experiment on 
+Jan 1st of year Y with yearly files and daily data to be interpolated, we do accept that the file 
+related to year Y-1 is not existing. The value of Jan 1st will be used as the missing one for 
+Dec 31st of year Y-1. If the file of year Y-1 exists, the code will read its last record. 
+Therefore, this file can contain only one record corresponding to Dec 31st, a useful feature for 
+user considering that it is too heavy to manipulate the complete file for year Y-1.
+
+
+% -------------------------------------------------------------------------------------------------------------
+% Interpolation on the Fly
+% -------------------------------------------------------------------------------------------------------------
+\subsection [Interpolation on-the-Fly] {Interpolation on-the-Fly}
+\label{SBC_iof}
+
+Interpolation on the Fly allows the user to supply input files required
+for the surface forcing on grids other than the model grid.
+To do this he or she must supply, in addition to the source data file,
+a file of weights to be used to interpolate from the data grid to the model grid.
+The original development of this code used the SCRIP package (freely available 
+\href{http://climate.lanl.gov/Software/SCRIP}{here} under a copyright agreement).
+In principle, any package can be used to generate the weights, but the
+variables in the input weights file must have the same names and meanings as
+assumed by the model.
+Two methods are currently available: bilinear and bicubic interpolation.
+
+\subsubsection{Bilinear Interpolation}
+\label{SBC_iof_bilinear}
+
+The input weights file in this case has two sets of variables: src01, src02,
+src03, src04 and wgt01, wgt02, wgt03, wgt04.
+The "src" variables correspond to the point in the input grid to which the weight
+"wgt" is to be applied. Each src value is an integer corresponding to the index of a 
+point in the input grid when written as a one dimensional array.  For example, for an input grid
+of size 5x10, point (3,2) is referenced as point 8, since (2-1)*5+3=8.
+There are four of each variable because bilinear interpolation uses the four points defining
+the grid box containing the point to be interpolated.
+All of these arrays are on the model grid, so that values src01(i,j) and
+wgt01(i,j) are used to generate a value for point (i,j) in the model.
+
+Symbolically, the algorithm used is:
+
+\begin{equation}
+f_{m}(i,j) = f_{m}(i,j) + \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}
+\end{equation}
+where function idx() transforms a one dimensional index src(k) into a two dimensional index,
+and wgt(1) corresponds to variable "wgt01" for example.
+
+\subsubsection{Bicubic Interpolation}
+\label{SBC_iof_bicubic}
+
+Again there are two sets of variables: "src" and "wgt".
+But in this case there are 16 of each.
+The symbolic algorithm used to calculate values on the model grid is now:
+
+\begin{equation*} \begin{split}
+f_{m}(i,j) =  f_{m}(i,j) +& \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}     
+              +   \sum_{k=5}^{8} {wgt(k)\left.\frac{\partial f}{\partial i}\right| _{idx(src(k))} }    \\
+              +& \sum_{k=9}^{12} {wgt(k)\left.\frac{\partial f}{\partial j}\right| _{idx(src(k))} }   
+              +   \sum_{k=13}^{16} {wgt(k)\left.\frac{\partial ^2 f}{\partial i \partial j}\right| _{idx(src(k))} }
+\end{split}
+\end{equation*}
+The gradients here are taken with respect to the horizontal indices and not distances since the spatial dependency has been absorbed into the weights.
+
+\subsubsection{Implementation}
+\label{SBC_iof_imp}
+
+To activate this option, a non-empty string should be supplied in the weights filename column 
+of the relevant namelist; if this is left as an empty string no action is taken.
+In the model, weights files are read in and stored in a structured type (WGT) in the fldread 
+module, as and when they are first required.
+This initialisation procedure determines whether the input data grid should be treated 
+as cyclical or not by inspecting a global attribute stored in the weights input file.
+This attribute must be called "ew\_wrap" and be of integer type.
+If it is negative, the input non-model grid is assumed not to be cyclic.
+If zero or greater, then the value represents the number of columns that overlap.
+$E.g.$ if the input grid has columns at longitudes 0, 1, 2, .... , 359, then ew\_wrap should be set to 0;
+if longitudes are 0.5, 2.5, .... , 358.5, 360.5, 362.5, ew\_wrap should be 2.
+If the model does not find attribute ew\_wrap, then a value of -999 is assumed.
+In this case the \rou{fld\_read} routine defaults ew\_wrap to value 0 and therefore the grid 
+is assumed to be cyclic with no overlapping columns.
+(In fact this only matters when bicubic interpolation is required.)
+Note that no testing is done to check the validity in the model, since there is no way 
+of knowing the name used for the longitude variable,
+so it is up to the user to make sure his or her data is correctly represented.
+
+Next the routine reads in the weights.
+Bicubic interpolation is assumed if it finds a variable with name "src05", otherwise 
+bilinear interpolation is used. The WGT structure includes dynamic arrays both for 
+the storage of the weights (on the model grid), and when required, for reading in 
+the variable to be interpolated (on the input data grid).
+The size of the input data array is determined by examining the values in the "src" 
+arrays to find the minimum and maximum i and j values required.
+Since bicubic interpolation requires the calculation of gradients at each point on the grid, 
+the corresponding arrays are dimensioned with a halo of width one grid point all the way around.
+When the array of points from the data file is adjacent to an edge of the data grid, 
+the halo is either a copy of the row/column next to it (non-cyclical case), or is a copy 
+of one from the first few columns on the opposite side of the grid (cyclical case).
+
+\subsubsection{Limitations}
+\label{SBC_iof_lim}
+
+\begin{enumerate}  
+\item  The case where input data grids are not logically rectangular has not been tested.
+\item  This code is not guaranteed to produce positive definite answers from positive definite inputs
+          when a bicubic interpolation method is used.
+\item  The cyclic condition is only applied on left and right columns, and not to top and bottom rows.
+\item  The gradients across the ends of a cyclical grid assume that the grid spacing between 
+          the two columns involved are consistent with the weights used.
+\item  Neither interpolation scheme is conservative. (There is a conservative scheme available 
+          in SCRIP, but this has not been implemented.)
+\end{enumerate}
+
+\subsubsection{Utilities}
+\label{SBC_iof_util}
+
+% to be completed
+A set of utilities to create a weights file for a rectilinear input grid is available 
+(see the directory NEMOGCM/TOOLS/WEIGHTS).
+
+
+% ================================================================
 % Analytical formulation (sbcana module) 
 % ================================================================
@@ -215,19 +470,5 @@
 read in the file, the time frequency at which it is given (in hours), and a logical 
 setting whether a time interpolation to the model time step is required 
-for this field). (fld\_i namelist structure).
-
-\textbf{Caution}: when the frequency is set to --12, the data are monthly 
-values. These are assumed to be climatological values, so time interpolation 
-between December the 15$^{th}$ and January the 15$^{th}$ is done using 
-records 12 and 1
-
-When higher frequency is set and time interpolation is demanded, the model 
-will try to read the last (first) record of previous (next) year in a file 
-having the same name but a suffix {\_}prev{\_}year ({\_}next{\_}year) being 
-added (e.g. "{\_}1989"). These files must only contain a single record. If they don't exist, 
-the model assumes that the last record of the previous year is equal to the first 
-record of the current year, and similarly, that the first record of the 
-next year is equal to the last record of the current year. This will cause 
-the forcing to remain constant over the first and last half fld\_frequ hours.
+for this field. See \S\ref{SBC_fldread} for a more detailed description of the parameters.
 
 Note that in general, a flux formulation is used in associated with a 
@@ -281,5 +522,5 @@
 \begin{table}[htbp]   \label{Tab_CORE}
 \begin{center}
-\begin{tabular}{|l|l|l|l|}
+\begin{tabular}{|l|c|c|c|}
 \hline
 Variable desciption					& Model variable	& Units	 & point \\		\hline
@@ -297,6 +538,8 @@
 %--------------------------------------------------------------------------------------------------------------
 
-Note that the air velocity is provided at a tracer ocean point, not at a velocity ocean point ($u$- and $v$-points). It is simpler and faster (less fields to be read), but it is not the recommended method when the ocean grid
-size is the same or larger than the one of the input atmospheric fields.
+Note that the air velocity is provided at a tracer ocean point, not at a velocity ocean 
+point ($u$- and $v$-points). It is simpler and faster (less fields to be read), 
+but it is not the recommended method when the ocean grid size is the same 
+or larger than the one of the input atmospheric fields.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -338,6 +581,5 @@
 As for the flux formulation, information about the input data required by the 
 model is provided in the namsbc\_blk\_core or namsbc\_blk\_clio 
-namelist (via the structure fld\_i). The first and last record assumption is also made 
-(see \S\ref{SBC_flx})
+namelist (see \S\ref{SBC_fldread}). 
 
 % ================================================================
@@ -399,5 +641,5 @@
 (see \mdl{dynspg} for the ocean). For sea-ice, the sea surface height, $\eta_m$, 
 which is provided to the sea ice model is set to $\eta - \eta_{ib}$ (see \mdl{sbcssr} module).
-$\eta_{ib}$ can be set in the output. This can simplify the altirmetry data and model comparison 
+$\eta_{ib}$ can be set in the output. This can simplify altimetry data and model comparison 
 as inverse barometer sea surface height is usually removed from these date prior to their distribution.
 
@@ -433,5 +675,6 @@
 River runoff generally enters the ocean at a nonzero depth rather than through the surface.
 Many models, however, have traditionally inserted river runoff to the top model cell.
-This was the case in \NEMO prior to the version 3.3, and was combined with an option to increase vertical mixing near the river mouth.
+This was the case in \NEMO prior to the version 3.3, and was combined with an option 
+to increase vertical mixing near the river mouth.
 
 However, with this method numerical and physical problems arise when the top grid cells are 
@@ -517,121 +760,4 @@
 
 }
-
-
-% ================================================================
-% Interpolation on the Fly
-% ================================================================
-
-\section [Interpolation on the Fly] {Interpolation on the Fly}
-\label{SBC_iof}
-
-Interpolation on the Fly allows the user to supply input files required
-for the surface forcing on grids other than the model grid.
-To do this he or she must supply, in addition to the source data file,
-a file of weights to be used to interpolate from the data grid to the model
-grid.
-The original development of this code used the SCRIP package (freely available 
-\href{http://climate.lanl.gov/Software/SCRIP}{here} under a copyright agreement).
-In principle, any package can be used to generate the weights, but the
-variables in the input weights file must have the same names and meanings as
-assumed by the model.
-Two methods are currently available: bilinear and bicubic interpolation.
-
-\subsection{Bilinear Interpolation}
-\label{SBC_iof_bilinear}
-
-The input weights file in this case has two sets of variables: src01, src02,
-src03, src04 and wgt01, wgt02, wgt03, wgt04.
-The "src" variables correspond to the point in the input grid to which the weight
-"wgt" is to be applied. Each src value is an integer corresponding to the index of a 
-point in the input grid when written as a one dimensional array.  For example, for an input grid
-of size 5x10, point (3,2) is referenced as point 8, since (2-1)*5+3=8.
-There are four of each variable because bilinear interpolation uses the four points defining
-the grid box containing the point to be interpolated.
-All of these arrays are on the model grid, so that values src01(i,j) and
-wgt01(i,j) are used to generate a value for point (i,j) in the model.
-
-Symbolically, the algorithm used is:
-
-\begin{equation}
-f_{m}(i,j) = f_{m}(i,j) + \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}
-\end{equation}
-where function idx() transforms a one dimensional index src(k) into a two dimensional index,
-and wgt(1) corresponds to variable "wgt01" for example.
-
-\subsection{Bicubic Interpolation}
-\label{SBC_iof_bicubic}
-
-Again there are two sets of variables: "src" and "wgt".
-But in this case there are 16 of each.
-The symbolic algorithm used to calculate values on the model grid is now:
-
-\begin{equation*} \begin{split}
-f_{m}(i,j) =  f_{m}(i,j) +& \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}     
-              +   \sum_{k=5}^{8} {wgt(k)\left.\frac{\partial f}{\partial i}\right| _{idx(src(k))} }    \\
-              +& \sum_{k=9}^{12} {wgt(k)\left.\frac{\partial f}{\partial j}\right| _{idx(src(k))} }   
-              +   \sum_{k=13}^{16} {wgt(k)\left.\frac{\partial ^2 f}{\partial i \partial j}\right| _{idx(src(k))} }
-\end{split}
-\end{equation*}
-The gradients here are taken with respect to the horizontal indices and not distances since the spatial dependency has been absorbed into the weights.
-
-\subsection{Implementation}
-\label{SBC_iof_imp}
-
-To activate this option, a non-empty string should be supplied in the weights filename column 
-of the relevant namelist; if this is left as an empty string no action is taken.
-In the model, weights files are read in and stored in a structured type (WGT) in the fldread 
-module, as and when they are first required.
-This initialisation procedure determines whether the input data grid should be treated 
-as cyclical or not by inspecting a global attribute stored in the weights input file.
-This attribute must be called "ew\_wrap" and be of integer type.
-If it is negative, the input non-model grid is assumed not to be cyclic.
-If zero or greater, then the value represents the number of columns that overlap.
-$E.g.$ if the input grid has columns at longitudes 0, 1, 2, .... , 359, then ew\_wrap should be set to 0;
-if longitudes are 0.5, 2.5, .... , 358.5, 360.5, 362.5, ew\_wrap should be 2.
-If the model does not find attribute ew\_wrap, then a value of -999 is assumed.
-In this case the \rou{fld\_read} routine defaults ew\_wrap to value 0 and therefore the grid 
-is assumed to be cyclic with no overlapping columns.
-(In fact this only matters when bicubic interpolation is required.)
-Note that no testing is done to check the validity in the model, since there is no way 
-of knowing the name used for the longitude variable,
-so it is up to the user to make sure his or her data is correctly represented.
-
-Next the routine reads in the weights.
-Bicubic interpolation is assumed if it finds a variable with name "src05", otherwise 
-bilinear interpolation is used. The WGT structure includes dynamic arrays both for 
-the storage of the weights (on the model grid), and when required, for reading in 
-the variable to be interpolated (on the input data grid).
-The size of the input data array is determined by examining the values in the "src" 
-arrays to find the minimum and maximum i and j values required.
-Since bicubic interpolation requires the calculation of gradients at each point on the grid, 
-the corresponding arrays are dimensioned with a halo of width one grid point all the way around.
-When the array of points from the data file is adjacent to an edge of the data grid, 
-the halo is either a copy of the row/column next to it (non-cyclical case), or is a copy 
-of one from the first few columns on the opposite side of the grid (cyclical case).
-
-\subsection{Limitations}
-\label{SBC_iof_lim}
-
-\begin{description}
-\item
-The case where input data grids are not logically rectangular has not been tested.
-\item
-This code is not guaranteed to produce positive definite answers from positive definite inputs.
-\item
-The cyclic condition is only applied on left and right columns, and not to top and bottom rows.
-\item
-The gradients across the ends of a cyclical grid assume that the grid spacing between the two columns involved are consistent with the weights used.
-\item
-Neither interpolation scheme is conservative.
-(There is a conservative scheme available in SCRIP, but this has not been implemented.)
-\end{description}
-
-\subsection{Utilities}
-\label{SBC_iof_util}
-
-% to be completed
-A set of utilities to create a weights file for a rectilinear input grid is available 
-(see the directory NEMOGCM/TOOLS/WEIGHTS).
 
 % ================================================================
