Index: /NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.bib
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.bib	(revision 10372)
+++ /NEMO/trunk/doc/latex/NEMO/main/NEMO_manual.bib	(revision 10373)
@@ -547,4 +547,15 @@
 PAGES = {1285--1297},
 DOI = {10.5194/gmd-8-1285-2015}
+}
+
+@ARTICLE{Breivik_al_JPO2014,
+   AUTHOR = {{\O}yvind Breivik and Peter A.E.M. Janssen and Jean-Raymond Bidlot},
+   YEAR = {2014},
+   TITLE = "{Approximate Stokes Drift Profiles in Deep Water}",
+   JOURNAL = {JPO},
+   VOLUME = {44},
+   NUMBER = {9},
+   DOI = {10.1175/JPO-D-14-0020.1.},
+   PAGES = {2433--2445, arXiv:1406.5039}
 }
 
@@ -1655,4 +1666,14 @@
 }
 
+@TECHREPORT{Janssen_al_TM13,
+ author = {P.A.E.M. Janssen and {\O}. Breivik and K. Mogensen 
+           and F. Vitart and  M. Balmaseda and J.B. Bidlot and 
+           S. Keeley and M. Leut-becher and L. Magnusson and F. Molteni},
+ title = {Air-Sea Interaction and Surface Waves},
+ year = {2013},
+ volume = {712},
+ institution = {ECMWF},
+}
+
 @ARTICLE{Jayne_St_Laurent_GRL01,
   author = {S.R. Jayne and L.C. {St. Laurent}},
@@ -2992,4 +3013,13 @@
   volume = {359},
   pages = {123--129}
+}
+
+@INCOLLECTION{Stokes_1847,
+  author    = {G.G. Stokes},
+  title     = {On the theory of oscillatory waves},
+  booktitle = {Transactions of the Cambridge Philosophy Society},
+  year      = {1847},
+  volume    = {8},
+  Pages     = {441--455}
 }
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 10372)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 10373)
@@ -32,12 +32,10 @@
 \end{itemize}
 
-Five different ways to provide the first six fields to the ocean are available which are controlled by
+Four different ways to provide the first six fields to the ocean are available which are controlled by
 namelist \ngn{namsbc} variables:
 an analytical formulation (\np{ln\_ana}\forcode{ = .true.}),
 a flux formulation (\np{ln\_flx}\forcode{ = .true.}),
 a bulk formulae formulation (CORE (\np{ln\_blk\_core}\forcode{ = .true.}),
-CLIO (\np{ln\_blk\_clio}\forcode{ = .true.}) or
-MFS \footnote { Note that MFS bulk formulae compute fluxes only for the ocean component}
-(\np{ln\_blk\_mfs}\forcode{ = .true.}) bulk formulae) and
+CLIO (\np{ln\_blk\_clio}\forcode{ = .true.}) bulk formulae) and
 a coupled or mixed forced/coupled formulation (exchanges with a atmospheric model via the OASIS coupler)
 (\np{ln\_cpl} or \np{ln\_mixcpl}\forcode{ = .true.}). 
@@ -76,7 +74,13 @@
   the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle
   (\np{ln\_dm2dc}\forcode{ = .true.});
-  and a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}\forcode{ = .true.}). 
+\item
+  a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}\forcode{ = .true.});
+\item
+  the Stokes drift rom an external wave model can be accounted (\np{ln\_sdw}\forcode{ = .true.}); 
+\item
+  the Stokes-Coriolis term can be included (\np{ln\_stcor}\forcode{ = .true.});
+\item
+  the surface stress felt by the ocean can be modified by surface waves (\np{ln\_tauwoc}\forcode{ = .true.}).
 \end{itemize}
-The latter option is possible only in case core or mfs bulk formulas are selected.
 
 In this chapter, we first discuss where the surface boundary condition appears in the model equations.
@@ -593,6 +597,6 @@
 % Bulk formulation
 % ================================================================
-\section[Bulk formulation {(\textit{sbcblk\{\_core,\_clio,\_mfs\}.F90})}]
-			{Bulk formulation {(\protect\mdl{sbcblk\_core}, \protect\mdl{sbcblk\_clio}, \protect\mdl{sbcblk\_mfs})}}
+\section[Bulk formulation {(\textit{sbcblk\{\_core,\_clio\}.F90})}]
+                        {Bulk formulation {(\protect\mdl{sbcblk\_core}, \protect\mdl{sbcblk\_clio})}}
 \label{sec:SBC_blk}
 
@@ -600,8 +604,8 @@
 
 The atmospheric fields used depend on the bulk formulae used.
-Three bulk formulations are available:
-the CORE, the CLIO and the MFS bulk formulea.
+Two bulk formulations are available:
+the CORE and the CLIO bulk formulea.
 The choice is made by setting to true one of the following namelist variable:
-\np{ln\_core} ; \np{ln\_clio} or  \np{ln\_mfs}.
+\np{ln\_core} or \np{ln\_clio}.
 
 Note:
@@ -712,41 +716,4 @@
 the namsbc\_blk\_core or namsbc\_blk\_clio namelist (see \autoref{subsec:SBC_fldread}). 
 
-% -------------------------------------------------------------------------------------------------------------
-%        MFS Bulk formulae
-% -------------------------------------------------------------------------------------------------------------
-\subsection{MFS formulea (\protect\mdl{sbcblk\_mfs}, \protect\np{ln\_mfs}\forcode{ = .true.})}
-\label{subsec:SBC_blk_mfs}
-%------------------------------------------namsbc_mfs----------------------------------------------------
-%
-%\nlst{namsbc_mfs}
-%----------------------------------------------------------------------------------------------------------
-
-The MFS (Mediterranean Forecasting System) bulk formulae have been developed by \citet{Castellari_al_JMS1998}. 
-They have been designed to handle the ECMWF operational data and are currently in use in
-the MFS operational system \citep{Tonani_al_OS08}, \citep{Oddo_al_OS09}.
-The wind stress computation uses a drag coefficient computed according to \citet{Hellerman_Rosenstein_JPO83}.
-The surface boundary condition for temperature involves the balance between
-surface solar radiation, net long-wave radiation, the latent and sensible heat fluxes.
-Solar radiation is dependent on cloud cover and is computed by means of an astronomical formula \citep{Reed_JPO77}.
-Albedo monthly values are from \citet{Payne_JAS72} as means of the values at $40^{o}N$ and $30^{o}N$ for
-the Atlantic Ocean (hence the same latitudinal band of the Mediterranean Sea).
-The net long-wave radiation flux \citep{Bignami_al_JGR95} is a function of
-air temperature, sea-surface temperature, cloud cover and relative humidity.
-Sensible heat and latent heat fluxes are computed by classical bulk formulae parameterised according to
-\citet{Kondo1975}.
-Details on the bulk formulae used can be found in \citet{Maggiore_al_PCE98} and \citet{Castellari_al_JMS1998}.
-
-Options are defined through the \ngn{namsbc\_mfs} namelist variables.
-The required 7 input fields must be provided on the model Grid-T and are:
-\begin{itemize}
-\item          Zonal Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windi})
-\item          Meridional Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windj})
-\item          Total Claud Cover (\%)  (\np{sn\_clc})
-\item          2m Air Temperature ($K$) (\np{sn\_tair})
-\item          2m Dew Point Temperature ($K$)  (\np{sn\_rhm})
-\item          Total Precipitation ${Kg} m^{-2} s^{-1}$ (\np{sn\_prec})
-\item          Mean Sea Level Pressure (${Pa}$) (\np{sn\_msl})
-\end{itemize}
-% -------------------------------------------------------------------------------------------------------------
 % ================================================================
 % Coupled formulation
@@ -1203,4 +1170,191 @@
 since its trajectory data may be spread across multiple files.
 
+% -------------------------------------------------------------------------------------------------------------
+%        Interactions with waves (sbcwave.F90, ln_wave)
+% -------------------------------------------------------------------------------------------------------------
+\section{Interactions with waves (\protect\mdl{sbcwave}, \protect\np{ln\_wave})}
+\label{sec:SBC_wave}
+%------------------------------------------namsbc_wave--------------------------------------------------------
+
+\nlst{namsbc_wave}
+%-------------------------------------------------------------------------------------------------------------
+
+Ocean waves represent the interface between the ocean and the atmosphere, so NEMO is extended to incorporate 
+physical processes related to ocean surface waves, namely the surface stress modified by growth and 
+dissipation of the oceanic wave field, the Stokes-Coriolis force and the Stokes drift impact on mass and 
+tracer advection; moreover the neutral surface drag coefficient from a wave model can be used to evaluate 
+the wind stress.
+
+Physical processes related to ocean surface waves can be accounted by setting the logical variable 
+\np{ln\_wave}\forcode{= .true.} in \ngn{namsbc} namelist. In addition, specific flags accounting for 
+different porcesses should be activated as explained in the following sections.
+
+Wave fields can be provided either in forced or coupled mode:
+\begin{description}
+\item[forced mode]: wave fields should be defined through the \ngn{namsbc\_wave} namelist 
+for external data names, locations, frequency, interpolation and all the miscellanous options allowed by 
+Input Data generic Interface (see \autoref{sec:SBC_input}). 
+\item[coupled mode]: NEMO and an external wave model can be coupled by setting \np{ln\_cpl} \forcode{= .true.} 
+in \ngn{namsbc} namelist and filling the \ngn{namsbc\_cpl} namelist.
+\end{description}
+
+
+% ================================================================
+% Neutral drag coefficient from wave model (ln_cdgw)
+
+% ================================================================
+\subsection{Neutral drag coefficient from wave model (\protect\np{ln\_cdgw})}
+\label{subsec:SBC_wave_cdgw}
+
+The neutral surface drag coefficient provided from an external data source ($i.e.$ a wave
+model), 
+can be used by setting the logical variable \np{ln\_cdgw} \forcode{= .true.} in \ngn{namsbc} namelist. 
+Then using the routine \rou{turb\_ncar} and starting from the neutral drag coefficent provided, 
+the drag coefficient is computed according to the stable/unstable conditions of the 
+air-sea interface following \citet{Large_Yeager_Rep04}. 
+
+
+% ================================================================
+% 3D Stokes Drift (ln_sdw, nn_sdrift)
+% ================================================================
+\subsection{3D Stokes Drift (\protect\np{ln\_sdw, nn\_sdrift})}
+\label{subsec:SBC_wave_sdw}
+
+The Stokes drift is a wave driven mechanism of mass and momentum transport \citep{Stokes_1847}. 
+It is defined as the difference between the average velocity of a fluid parcel (Lagrangian velocity) 
+and the current measured at a fixed point (Eulerian velocity). 
+As waves travel, the water particles that make up the waves travel in orbital motions but 
+without a closed path. Their movement is enhanced at the top of the orbit and slowed slightly 
+at the bottom so the result is a net forward motion of water particles, referred to as the Stokes drift. 
+An accurate evaluation of the Stokes drift and the inclusion of related processes may lead to improved 
+representation of surface physics in ocean general circulation models.
+The Stokes drift velocity $\mathbf{U}_{st}$ in deep water can be computed from the wave spectrum and may be written as: 
+
+\begin{equation} \label{eq:sbc_wave_sdw}
+\mathbf{U}_{st} = \frac{16{\pi^3}} {g} 
+                \int_0^\infty \int_{-\pi}^{\pi} (cos{\theta},sin{\theta}) {f^3}
+                \mathrm{S}(f,\theta) \mathrm{e}^{2kz}\,\mathrm{d}\theta {d}f
+\end{equation}
+
+where: ${\theta}$ is the wave direction, $f$ is the wave intrinsic frequency, 
+$\mathrm{S}($f$,\theta)$ is the 2D frequency-direction spectrum, 
+$k$ is the mean wavenumber defined as: 
+$k=\frac{2\pi}{\lambda}$ (being $\lambda$ the wavelength). \\
+
+In order to evaluate the Stokes drift in a realistic ocean wave field the wave spectral shape is required 
+and its computation quickly becomes expensive as the 2D spectrum must be integrated for each vertical level. 
+To simplify, it is customary to use approximations to the full Stokes profile.
+Three possible parameterizations for the calculation for the approximate Stokes drift velocity profile 
+are included in the code through the \np{nn\_sdrift} parameter once provided the surface Stokes drift 
+$\mathbf{U}_{st |_{z=0}}$ which is evaluated by an external wave model that accurately reproduces the wave spectra 
+and makes possible the estimation of the surface Stokes drift for random directional waves in 
+realistic wave conditions:
+
+\begin{description}
+\item[\np{nn\_sdrift} = 0]: exponential integral profile parameterization proposed by 
+\citet{Breivik_al_JPO2014}:
+
+\begin{equation} \label{eq:sbc_wave_sdw_0a}
+\mathbf{U}_{st} \cong \mathbf{U}_{st |_{z=0}} \frac{\mathrm{e}^{-2k_ez}} {1-8k_ez} 
+\end{equation}
+
+where $k_e$ is the effective wave number which depends on the Stokes transport $T_{st}$ defined as follows:
+
+\begin{equation} \label{eq:sbc_wave_sdw_0b}
+k_e = \frac{|\mathbf{U}_{\left.st\right|_{z=0}}|} {|T_{st}|} 
+\quad \text{and }\
+T_{st} = \frac{1}{16} \bar{\omega} H_s^2 
+\end{equation}
+
+where $H_s$ is the significant wave height and $\omega$ is the wave frequency.
+
+\item[\np{nn\_sdrift} = 1]: velocity profile based on the Phillips spectrum which is considered to be a 
+reasonable estimate of the part of the spectrum most contributing to the Stokes drift velocity near the surface
+\citep{Breivik_al_OM2016}:
+
+\begin{equation} \label{eq:sbc_wave_sdw_1}
+\mathbf{U}_{st} \cong \mathbf{U}_{st |_{z=0}} \Big[exp(2k_pz)-\beta \sqrt{-2 \pi k_pz} 
+\textit{ erf } \Big(\sqrt{-2 k_pz}\Big)\Big]
+\end{equation}
+
+where $erf$ is the complementary error function and $k_p$ is the peak wavenumber.
+
+\item[\np{nn\_sdrift} = 2]: velocity profile based on the Phillips spectrum as for \np{nn\_sdrift} = 1 
+but using the wave frequency from a wave model.
+
+\end{description}
+
+The Stokes drift enters the wave-averaged momentum equation, as well as the tracer advection equations 
+and its effect on the evolution of the sea-surface height ${\eta}$ is considered as follows: 
+
+\begin{equation} \label{eq:sbc_wave_eta_sdw}
+\frac{\partial{\eta}}{\partial{t}} = 
+-\nabla_h \int_{-H}^{\eta} (\mathbf{U} + \mathbf{U}_{st}) dz 
+\end{equation}
+
+The tracer advection equation is also modified in order for Eulerian ocean models to properly account 
+for unresolved wave effect. The divergence of the wave tracer flux equals the mean tracer advection 
+that is induced by the three-dimensional Stokes velocity. 
+The advective equation for a tracer $c$ combining the effects of the mean current and sea surface waves 
+can be formulated as follows: 
+
+\begin{equation} \label{eq:sbc_wave_tra_sdw}
+\frac{\partial{c}}{\partial{t}} = 
+- (\mathbf{U} + \mathbf{U}_{st}) \cdot \nabla{c}
+\end{equation}
+
+
+% ================================================================
+% Stokes-Coriolis term (ln_stcor)
+% ================================================================
+\subsection{Stokes-Coriolis term (\protect\np{ln\_stcor})}
+\label{subsec:SBC_wave_stcor}
+
+In a rotating ocean, waves exert a wave-induced stress on the mean ocean circulation which results 
+in a force equal to $\mathbf{U}_{st}$×$f$, where $f$ is the Coriolis parameter. 
+This additional force may have impact on the Ekman turning of the surface current. 
+In order to include this term, once evaluated the Stokes drift (using one of the 3 possible 
+approximations described in \autoref{subsec:SBC_wave_sdw}), 
+\np{ln\_stcor}\forcode{ = .true.} has to be set.
+
+
+% ================================================================
+% Waves modified stress (ln_tauwoc, ln_tauw)
+% ================================================================
+\subsection{Wave modified sress (\protect\np{ln\_tauwoc, ln\_tauw})} 
+\label{subsec:SBC_wave_tauw}
+
+The surface stress felt by the ocean is the atmospheric stress minus the net stress going 
+into the waves \citep{Janssen_al_TM13}. Therefore, when waves are growing, momentum and energy is spent and is not 
+available for forcing the mean circulation, while in the opposite case of a decaying sea 
+state more momentum is available for forcing the ocean. 
+Only when the sea state is in equilibrium the ocean is forced by the atmospheric stress, 
+but in practice an equilibrium sea state is a fairly rare event. 
+So the atmospheric stress felt by the ocean circulation $\tau_{oc,a}$ can be expressed as: 
+
+\begin{equation} \label{eq:sbc_wave_tauoc}
+\tau_{oc,a} = \tau_a - \tau_w
+\end{equation}
+
+where $\tau_a$ is the atmospheric surface stress;
+$\tau_w$ is the atmospheric stress going into the waves defined as:
+
+\begin{equation} \label{eq:sbc_wave_tauw}
+\tau_w = \rho g \int {\frac{dk}{c_p} (S_{in}+S_{nl}+S_{diss})}
+\end{equation}
+
+where: $c_p$ is the phase speed of the gravity waves,
+$S_{in}$, $S_{nl}$ and $S_{diss}$ are three source terms that represent 
+the physics of ocean waves. The first one, $S_{in}$, describes the generation 
+of ocean waves by wind and therefore represents the momentum and energy transfer 
+from air to ocean waves; the second term $S_{nl}$ denotes 
+the nonlinear transfer by resonant four-wave interactions; while the third term $S_{diss}$ 
+describes the dissipation of waves by processes such as white-capping, large scale breaking 
+eddy-induced damping.
+
+The wave stress derived from an external wave model can be provided either through the normalized 
+wave stress into the ocean by setting \np{ln\_tauwoc}\forcode{ = .true.}, or through the zonal and 
+meridional stress components by setting \np{ln\_tauw}\forcode{ = .true.}.
+
 
 % ================================================================
@@ -1421,24 +1575,4 @@
 \end{description}
 
-% -------------------------------------------------------------------------------------------------------------
-%        Neutral Drag Coefficient from external wave model
-% -------------------------------------------------------------------------------------------------------------
-\subsection[Neutral drag coeff. from external wave model (\protect\mdl{sbcwave})]
-				{Neutral drag coefficient from external wave model (\protect\mdl{sbcwave})}
-\label{subsec:SBC_wave}
-%------------------------------------------namwave----------------------------------------------------
-
-\nlst{namsbc_wave}
-%-------------------------------------------------------------------------------------------------------------
-
-In order to read a neutral drag coefficient, from an external data source ($i.e.$ a wave model),
-the logical variable \np{ln\_cdgw} in \ngn{namsbc} namelist must be set to \forcode{.true.}.
-The \mdl{sbcwave} module containing the routine \np{sbc\_wave} reads the namelist \ngn{namsbc\_wave}
-(for external data names, locations, frequency, interpolation and all the miscellanous options allowed by
-Input Data generic Interface see \autoref{sec:SBC_input}) and
-a 2D field of neutral drag coefficient.
-Then using the routine TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided, 
-the drag coefficient is computed according to stable/unstable conditions of the air-sea interface following
-\citet{Large_Yeager_Rep04}.
 
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/introduction.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/introduction.tex	(revision 10372)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/introduction.tex	(revision 10373)
@@ -311,4 +311,5 @@
 \item new definition of configurations ;
 \item bulk formulation ;
+\item NEMO-wave large scale interactions ;
 \item ... ; 
 \end{enumerate}
